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ABSTRACT 

We develop an algorithm for estimating parameters of a distribution sampled with contamination. 
We employ a statistical technique known as "expectation maximization" (EM). Given models for both 
member and contaminant populations, the EM algorithm iteratively evaluates the membership prob- 
ability of each discrete data point, then uses those probabilities to update parameter estimates for 
member and contaminant distributions. The EM approach has wide applicability to the analysis of 
astronomical data. Here we tailor an EM algorithm to operate on spectroscopic samples obtained with 
the Michigan-MIKE Fiber System (MMFS) as part of our Magellan survey of stellar radial velocities in 
nearby dwarf spheroidal (dSph) galaxies. These samples, to be presented in a companion paper, con- 
tain discrete measurements of line-of-sight velocity, projected position, and pseudo-equivalent width 
of the Mg-triplet feature, for ~ 1000 — 2500 stars per dSph, including some fraction of contamina- 
tion by foreground Milky Way stars. The EM algorithm uses all of the available data to quantify 
dSph and contaminant distributions. For distributions (e.g., velocity and Mg-index of dSph stars) 
assumed to be Gaussian, the EM algorithm returns maximum-likelihood estimates of the mean and 
variance, as well as the probability that each star is a dSph member. These probabilities can serve as 
weights in subsequent analyses. Applied to our MMFS data, the EM algorithm identifies more than 
5000 stars as probable dSph members. We test the performance of the EM algorithm on simulated 
data sets that represent a range of sample size, level of contamination, and amount of overlap between 
dSph and contaminant velocity distributions. The simulations establish that for samples ranging from 
large {N ~ 3000, characteristic of the MMFS samples) to small {N ~ 30, resembling new samples 
for extremely faint dSphs), the EM algorithm distinguishes members from contaminants and returns 
accurate parameter estimates much more reliably than conventional methods of contaminant removal 
(e.g., Sigma clipping). 

Subject headings: galaxies: dwarf — galaxies: kinematics and dynamics — (galaxies:) Local Group 
— galaxies: individual (Carina, Fornax, Sculptor, Sextans) — techniques: radial 
velocities 



1. INTRODUCTION 

Most astronomical data sets are polluted to some 
extent by foreground/background objects ("contami- 
nants" ) that can be difficult to distinguish from objects 
of interest ("members"). Contaminants may have the 
I same apparent magnitudes, colors, and even velocities 
as members, and so satisfy many of the criteria used 
to identify members prior to and/or after observation. 
Obviously, all else being equal, analyses in which those 
contaminants are properly identified are superior to anal- 
yses in which they are not. Typically, one attempts to 
remove contaminants prior to analysis by assuming some 
reasonable distribution (e.g., Gaussian) for the member 
population, and then rejecting outliers iteratively until 
the number of outliers and the parameters (e.g., mean 
and variance) of the assumed distribution stabilize. 

This conventional approach is problematic in that the 
results depend on the arbitrary definition of some thresh- 
old (e.g., 3(t) separating members from outliers. Further, 
as they do not incorporate any explicit consideration of 
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the contaminant distribution, conventional methods pro- 
vide no means of identifying likely contaminants that lie 
near the center of the member distribution. The resulting 
samples are thus likely to retain contaminants that will 
then receive equal weight as true members in subsequent 
analyses. 

In this paper we promote a different approach to the 
problem of sample contamination, drawing upon a tech- 
nique known in the s tatistics literature as "expectation 
maximization" (EM: 'Demuster et al.lfl977l ). We intro- 
duce an iterative EM algorithm for estimating the mean 
and variance of a distribution sampled with contamina- 
tion. The EM method differs from a conventional, sigma- 
clipping method in two key respects. First, whereas the 
conventional method answers the question of member- 
ship with a yes or no, the EM method yields a prob- 
ability of membership. This probability serves as a 
weight during subsequent analysis. No data points are 
discarded; rather, likely contaminants receive appropri- 
ately little weight. Second, in evaluating membership 
probability, the EM method explicitly considers the dis- 
tributions of both members and contaminants. Thus 
likely contaminants that happen to lie near the center 
of the member distribution can be identified as such 
and weighted accordingly. EM therefore allows a full 
maximum- likelihood treatment of the entire data set. 
The EM approach has potentially wide applicability to 
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many types of data. In previous work, ISen et al] ()2007f ) 
briefly discuss the application of an EM algorithm to 
dSph data as a means for separating members from con- 
taminants. Here we exapand upon this work and develop 
an EM algorithm in order to estimate, given contami- 
nated spectroscopic samples, mean velocities and veloc- 
ity dispersions of press ure-supported galaxies. 

In companion papers (jWalker et al.ll2007ll2008l Papers 
I and II, respectively), we present a spectroscopic sample 
from our Michigan/MIKE Fiber System (MMFS) sur- 
vey of stellar radial velocities in nearby dwarf spheroidal 
(dSph) galaxies. For each star we measure the line-of- 
sight velocity and the pseudo-equivalent width of the 
Mg-triplet absorption feature. As of 2008 August, the 
MMFS data set contains measurements of 7103 red gi- 
ant candidates in the dSphs Carina, Fornax, Sculptor 
and Sextans, each a satellite of the Milky Way. Along 
the line of sight to each of the targeted dSphs, the Milky 
Way contributes stars with magnitudes and colors sat- 
isfying our CMD-based target selection criteria (Figure 
1 of Paper I). Therefore we expect the MMFS sample 
to carry a degree of contamination that varies with the 
density contrast between dSph and foreground. 

Figures [1] and [5] plot for all MMFS stars the measured 
velocity, V, magnesium index SMg, and projected ra- 
dius R (angular distance from the dSph center). Viewing 
these scatter plots, one can distinguish loci of foreground 
stars from those of the Carina, Sculptor and Sextans pop- 
ulations. The members of these three dSphs cluster into 
relatively narrow velocity distributions and have system- 
atically weaker EMg than do the foreground stars. The 
latter effect owes to the fact that foreground stars suffi- 
ciently faint to be included in our CMD selection regions 
are likely to be K-dwarfs in the MW disk, which exhibit 
strong magnesium absorption due to the fact that their 
higher surface gravities and highe r metal ab undances 
increase their atmospheric opacity (jCavrel et aLlll99ll ). 
Due to Fornax's relatively high surface brightness, our 
Fornax targets are more likely to be bona fide members, 
and the Fornax data do not show an obviously distinct 
distribution of foreground stars. However, given Fornax's 
relatively broad distribution of SMg and small systemic 
velocity of ~ 55 km s^"'^ along the line of sight, its ve- 
locity distribution overlaps with that expected for the 
foreground distribution. Therefore we must not consider 
the Fornax sample to be free of contamination, especially 
in the low- velocity tail of the distribution (see Figure [T]). 

EM provides a reliable and efficient means of character- 
izing the distributions of both member and contaminant 
populations. In what follows we provide a general de- 
scription of the EM approach, which has many potential 
applications in astronomy. We then develop an EM algo- 
rithm for our specific purpose of measuring dSph mean 
velocities and velocity dispersions from the MMFS data. 
Our EM algorithm incorporates stellar velocity, magne- 
sium index and projected position to generate maximum- 
likelihood estimates of the desired parameters. The al- 
gorithm also evaluates for each star the probability of 
dSph membership. Adding these probabilities, we find 
that the MMFS sample contains more than 5000 dSph 
members. Finally, we generate artificial dSph-like sam- 
ples that allow us to examine the performance of the 
EM algorithm given a variety of sample sizes, levels of 
contamination, and amounts of overlap between member 



and contaminant velocity distributions. We demonstrate 
that the EM algorithm consistently recovers the correct 
distribution parameters in all but the smallest and most 
severely contaminated samples, and dramatically outper- 
forms conventional methods for outlier rejection. 

2. EXPECTATION-MAXIMIZATION 

Consider a population for which some random variable 
X follows a probability distribution Pmem (x) that is char- 
acterized by parameter set Cmem- For example, if the dis- 



tribution is Gaussian, Pmem{x) = (27r(T^) ^/^exp[— (x — 

(X))2/(2cr2)], and Cmem = {(-^)mem, CT^em} COUSistS of 

the mean and variance. Typically one samples the dis- 
tribution with the goal of estimating Cmem- If contam- 
ination is present, then some fraction of the sample is 
drawn from a second, "non-member" population that 
has the probability distribution Pnon{x), with parame- 
ter set Cnon- Let p{a) represent the unconditional (i.e., a 
priori)membei fraction which, unless constant over the 
sample, depends on some random variable A that is inde- 
pendent of X. For example, in the case of dSph data, we 
find that a star's a priori probability of membership cor- 
relates with its projected distance from the dSph center, 
such that the member fraction, p(a), is highest at the 
center of the galaxy and decreases toward its sparsely 
populated outer regions. 

Suppose the variable M indicates membership and 
takes one of two possible values: M = 1 indicates ob- 
servation of a member, M = indicates observation of a 
nonmember. Given a data set {Xi, Mi}^^ for which all 
Ali are known, one can evaluate the likelihood, 
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or more conveniently, the log-likelihood given by 
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However, the values Mj, which can be cither exactly 
or exactly 1, usually are unknown, in which case the 
likelihood cannot be evaluated. Using the EM technique, 
one evaluates instead the expected log-likelihood in terms 
of the expected values of M . Prior to observation, M is 
a Bernoulli random variable with "probability of sucess" 
given by the a priori membership fraction: P[M = 1) = 
p{a). After observation, the expected value of M is the 
probability that M = 1, subject to the data as well as 
the prior constraints specified by p(a). Denoting this 
probability Pm, we have 
Pm. =P(A/, = l|X„aO 

i{Xi)p{ai) 



_ Pn 

pnon {Xi)[l-p{a,)]- 

The expected log-likelihood, given data S = {Xi\f^ 
then 

N 
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Fig. 1. — Projected distance from the dSph center (top panels) and magnesium index (bottom panels) versus velocity, Carina and Fornax dSphs. 
Marker color indicates probability that the star belongs to the dSph population, evaluated using the EM algorithm. Black (red; magenta; green; 
cyan; blue) markers signify Pm ^ 0.01 {Pm > 0.01; > 0.50; > 0.68; > 0.95; > 0.99). In each main panel a pair of vertical, dotted lines encloses 
stars that would pass conventional membership tests based on an iterative, 3fT velocity threshold. Sub-panels at bottom give the observed velocity 



distribution. The red curve in the bottom is the expected velocity distribution of foreground contaminants, 
model. The predicted Besangon distribution along the line of sight to each dSph is normalized according 
Nm/N. 



derived from the Besangon Milky Way 
to the estimated membership fraction 
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In the EM approach, maximum-hkehhood estimates of 
the parameter sets Cmem and Cnon, as well as any un- 
known parameters in p{a), take the values that maximize 

£/(lni(^mem, Cnon)\S). 

The iterative EM algorithm derives its name from its 
two fundamental steps. In the "expectation" step, the 
expected value of Mi is identified with the probability 
Pm; and evaluated for all i according to Equation [3l 
In the "maximization" step, the distribution parameters 
are estimated by maximizing the expected log-likelihood 
specified by Equation [3) The parameter estimates are 
then used to update the values P^f. in the next expec- 
tation step, followed by the updating of parameter esti- 
mates in the subsequent maximization step, and so forth 



until convergence occurs, 
parameter sets Cmem, Crw 



In addition to estimates of the 
,„ (and any free parameters in 
the function p{a)), the algorithm provides useful esti- 
mates of Pm;, the probability of membership, for all i. 



3. EM APPLIED TO DSPH DATA 

We now develop an EM algorithm for estimating pa- 
rameters of the distributions sampled by the MMFS 
data. In this application we evaluate the membership 
of stars using MMFS data that sample two independent 
distributions, velocity and magnesium index. Also, in- 
spection of the top panels in Figures [1] and reveals 
that the fraction of contamination tends to increase as a 
function of distance from the dSph center. Projected ra- 
dius therefore serves as the independent variable in the 
unconditional probability function p{a). In a previous 
application of the EM algorithm to dSph data, ISen et al] 
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Fig. 2. — Same as Figuicf^ but for the Sculptor and Sextans dSphs. 



assume p{a) decreases exponentially, consistent 
with the measured surface br ightness profiles of dSphs 
(jlrwin fc Hatzidimitrioul [19951 ). Here, in order to avoid 
bias due to selection effects, we choose not to specify a 
particular functional form for p{a); instead, after each 
maximization step we estimate p{a) via monotonic re- 
gression (see Appendix |BJ . In the end we obtain esti- 
mates of the means and variances of the velocity and 
EMg distributions for member and nonmember popula- 
tions, the membership probability of each star, and the 
member fraction as a function of projected radius. 

We assume that for dSph members the joint distri- 
bution of velocity, V, and magnesium strength, hence- 
forth denoted W, is a bivariate Gaussian distribution 



with means {V)mem and {W)r 



variances cr 



and 



Wo ,mem 



, and covariance avo 



Wo ,rnem 



Vq .mem 

{{V - 



{V)mem){W — {W)mem)) = (the distributions of V and 
W are independent) . If we were to sample only the dSph 
population (i.e., with no contamination), the probability 
density for observation of a star with velocity Vi and line 



strength Wi would be 



exp 



Pmemiy^, Wi) = 



^0 . 



■{V)r 



[Wi- 



'Wo,.! 



+ <^vM^Wo,mem+^w) 

(5) 

where and awi represent measurement errors. For 
the ~ 5% of stars in the MMFS sample that lack an ac- 
ceptable measurement of magnesium strength we replace 
Equation [5] with the univariate Gaussian probability 



exp 



Vq ,mem 



(6) 



A virtue of the EM method lies in its ability to incorpo- 
rate and to evaluate the contaminant as well as member 
distributions. The efficacy of the algorithm therefore de- 
pends on the suitability of the foreground model used 
to characterize Pnon{Vi,Wi). Since the distribution of 
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the MW's disk rotational velocities is non-Gaussian, we 
do not model the contaminant velocity distribution as a 
simple Gaussian. Instead we adopt the distribution of 
foreground velocities specified by the B esangon numeri- 
cal model of the MW ( Robin et a l.'2003')^ From the Be- 
sangon model, which includes both disk and halo compo- 
nents, we obtain simulated velocity and V, I photometric 
data along the line of sight to each dSph. After remov- 
ing simulated stars that fall outside our CMD selection 
regions (see Figure 1 of Paper I) , we use a Gaussian ker- 
nel to estimate the marginal distribution of contaminant 
velocities from Ni,es artificial data points: 



Pnon{v) = Pbesiv) 



Nbe. 



X exp 



1 jVbes 

2 a 



(7) 



We adopt^ "'Vtea = 2 km s^^, similar to the median 
MMFS velocity error. Note that this smooth foreground 
model does not account for known and/or potential sub- 
structu res (e.g. , debris from accretion events as discussed 
e.g., bv lBell et al. 2008) in the MW halo; for some appli- 
cations it may become necessary to build such features 
into the foreground model. 

We approximate the marginal distribution of fore- 
ground magnesium strength as Gaussian and indepen- 
dent of velocity, with mean {W)non and variance a"^^ 
Notice that this assumption ignores the fact that the 
MW foreground has distinct components — thin and thick 
disk, halo — with different Mg distributions. In some sit- 
uations, particularly when samples are small, it may 
become necessary to adopt a more realistic foreground 
model, but for our present purposes the single-Gaussian 
model suffices. Under this simple model, if we were to 
sample only from the population of contaminating MW 
stars satisfying our target-selection criteria, the proba- 
bility density for observation of a star with velocity Vi 
and line strength Wi would then be 



PbesiVi) exp 



PnoniViyWi) = 
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In reality we sample from the union of dSph and con- 
taminant populations. Assuming the surface brightness 
of nonmembers remains approximately constant over the 
face of a given dSph, the a priori probability of observ- 
ing a dSph member decreases in proportion to the surface 
density of actual members meeting our target selection 
criteria. We therefore choose the function p{a) to repre- 
sent the fraction of selected targets at elliptical radius^ 

^ available at http://bison.obs-besancon.fr/modole/ 
® We experimented with different choices for the smoothing 
bandwidth, including o"v^^, = 6 km as well as a variable band- 
width with each artificial error drawn randomly from the set of 
MMFS velocity errors. The EM algorithm converged on the same 
(to within 0.001 km s~^) parameter estimates and membership 
probabilities in all cases. 

A star's "elliptical radius" refers to the semi-major axis of 
the isophotal ellipse that passes through the star's position. We 
make the approximation that a dSph's isophotal ellipses all have 
common centers an d orientation, and we adopt e llipticities and 
position angles from llrwin &: Hatzidimitrioul I I1995I) . 



a that are actually dSph members. One might reason- 
ably expect p{a) to decrease exponentially, since expo- 
nential profiles generally provide adequate fits to dSph 
surface brightness data (jlrwin fc Hatzidimitrioul Il995l. 
However, in order to allow for sample bias introduced 
by our target selection procedure, we adopt the less re- 
strictive assumption that p{a) is merely a non-increasing 
function. 

Thus in our formulation, the EM analysis of a dSph for 
which we have data set S = {Vi,Wi,ai}^i involves the 
estimation of a set of 

{(^)mem; ^Vo,mem^ {W)memi ^Wo,7nem-! 

as well as the values of p{ai) for all i. Adopting the 

notation that X is the estimate of X, we let C^"> 
denote the set of parameter estimates obtained in the 
n*'* iteration of the algorithm. In the "expectation" 
step of the next iteration, we use these estimates to 



obtain pmJ^}, pk^tP , and then 
equations El [H and [3] we obtain 
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(11) 

In the "maximization" step of iteration n -I- 1 we max- 
imize the expected log-likelihood (Equation 2]). Because 
the form of p(a) is unknown, our maximization step has 
two parts. In the first part we estimate the parameter set 
C. By setting equal to zero the partial derivatives of the 
expected log-likelihood with respect to each parameter, 
we obtain six equations: 
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= 0: 



= 0: 



These equations do not have analytic solutions, but we 
obtain consistent estimates of the means and variances 
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using the estimates from the previous iteration. For ex- 
ample, 
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Iterative solutions for the means and variances of the 
magnesium distributions take the same form (with (1 — 



pj;"-'^) replacing P^j ' for estimates pertaining to the con- 
taminant distribution). We calculate the error associated 
with each parameter estimate by propagating the mea- 
surement error as well as the parameter errors from the 
previous iteration. Formulae for computing the sizes of 
1(7 errorbars are derived in Appendix lAl 

Recall that we assume about p{a) only that it is non- 
increasing. Therefore in the second part of the maxi- 
mization step, we maximize the expected log-likelihood 
(Equ ation with resp e ct to all non-increasing functions 
p{a). [Robertson et al.l (|1988 ) show that the solution is 
the isotonic estimator of the form 
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{n+l} 



(a») 



mm 

l<M<j 



max 

i<v<N V ■ 
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(15) 



where the notation maxa<z<t, f(z) specifies the maxi- 
mum value of f(z) in the interval a < z < b, and 
miTic<x<d g{x) specifies the minimum value of g{x) in the 
interval c < x < d. It is computationally expensive to 
perform an exhaustive search for the minimum among 
the maxima at each data point; however, the mono- 
tonic regression is performed effi ciently using the "Pool- 
Adjacent- Violators" algorithm (jGrotzinger fc Witzgalll 
[l98l . which we describe in Appendix iBl 

4. PROCEDURE 
In our implementation of the EM algorithm, we ini- 
tialize p{°>(ai) = 0.5 and P|°/ = 0.5 for aU i. We then 
initialize the variances a"^^ 

.mem ' ,mem 

and cr^^ 

and evaluate Equations [13] and [14] to obtain intial esti- 
mates of all parameters As Figures [3] - [6] show, 
the final estimates are insensitive to the values used to 
initialize the variances. 

After initialization, the algorithm begins with the ex- 
pectation step. We use Equations [9l - fTT] and the initial 

estimates in C'^"'' to calculate pml„i(Fi, Wi), pk])niVi,Wi), 
and Pj^j^ for all i. Next, in the maximization step, we 
use Equations [T2] and [TJ] to update the parameter es- 
timates in ^^^^ ■ After maximization, we use monotonic 
regression (see Appendix [B|) to evaluate Equation[T5l the 
regression provides estimates p^^^{ai) for all i. We then 
proceed to the next iteration, which begins again with 
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Fig. 3.— 

obtained in the first 15 iterations of the EM algorithm, applied here 
to tlie MMFS sample for the Carina dSph. In each panel, solid curves 
represent the member population and dotted curves represent the 
nonmember population. We have run the algorithm three times, 

using initial values {(TY^^rnerm ^WQ,merm ^WQ,non} = {5,0.1,0.1}, 

{50, 0.5, 0.5}, and {100, 1.0, 1.0} in units of {kms-\ A, A}. The final 
estimates are insensitive to the initialization, even for the large range 
of values considered. For illustrative purposes only, the velocity mean 
and dispersion shown for nonmembers are calculated as in Equations 
1131 the nonmember velocity distribution is non-Gaussian, and during 
the EM algorithm Pnoni^) is evaluated according to Equation [7] fsee 
Section [31 ■ 

the expectation step. Parameter estimates typically sta- 
bilize after ^ 15 — 20 iterations. Because the algorithm 
runs quickly, we iterate 50 times. 

Figures [3] - [S] depict the evolution of the parameters 
in Q"^^ for each dSph and for several choices of variance 
initialization. The estimates in Q"^^ obtained after 15-20 
iterations do not depend on the initialization. 

5. RESULTS 

Table [1] summarizes results from our application of the 
EM algorithm to the MMFS data for each dSph. Col- 
umn 2 lists the number of observed stars, N , and col- 
umn 3 lists the number of members, defined by Nm = 
'l2iLi ■ Columns 4-9 then list the parameter esti- 
mates in C. 

The EM algorithm effectively separates likely mem- 
bers from nonmembers in the MMFS data. Applying 
the algorithm to our entire MMFS sample, we find that 
4946 stars have Pm > 0.9, 1918 stars have Pm < 0.1, 
and 239 stars have Pm between these limits. Just 97 
stars have what we might call "ambiguous" membership, 
0.25 < Pm < 0.75. Due to the bimodal distribution of 
Pm values, Nm typically exceeds by less than 1% the 
number of stars having Pm > 0.9. Adding all proba- 
bilities, the entire MMFS data set contains Nm = 5063 
dSph member stars. 

Marker colors in Figures [T] and [2] indicate the values 
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Fig. 4. — Same as Figurc[3]but for Fornax data. 



Fig. 



Same as Figure [3] but for Sextans data 
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Fig. 5. — Same as Figure[3]but for Sculptor data. 

of the membership probabilities, Pm, returned by the 
EM algorithm. Black points signify the most likely non- 
members, with Pm < 0.01. For the remaining stars, 
bluer colors signify more likely membership: red, ma- 
genta, green, cyan, and blue points correspond to stars 
having Pm exceeding 0.01, 0.50, 0.68, 0.95, and 0.99, 
respectively. Pairs of dotted, vertical lines enclose stars 
that would be accepted as members using a conventional, 
3(7 velocity threshold (Section l7.ip . The nonmember ve- 
locity distributions predicted by the Bcsangon models are 
normalized for each dSph according to the final member- 



ship fraction and plotted in red in the bottom subpanels 
in Figures [1] and [51 We find that the velocity distribu- 
tions of the observed nonmembers generally agree with 
the Besangon predictions. 

Upper sub-panels in Figure [7| plot the measured veloc- 
ities as a function of elliptical radius. Bottom sub-panels 
plot the monotonic regression estimate p{a), the mem- 
ber fraction as a function of elliptical radius. We note 
that p{a) need not have a steep slope. Our simple as- 
sumption that p{a) is non-increasing allows the data to 
determine the dependence of membership on position. 
Figure [5] maps sky positions of the observed stars, again 
using marker color to indicate membership. 



6. 



PERFORMANCE OF THE EM ALGORITHM USING 
SIMULATED DATA 



In order to evaluate the performance of the EM al- 
gorithm under a variety of conditions, we test it on 
simulated data sets that differ in sample size, level of 
contamination, and amount of velocity overlap between 
member and contaminant distributions. In all simula- 
tions we draw a fraction fmem of N artificial data points 
from a dSph member population that we assume to have 
a Gaussian velocity distribution with dispersion 10 km 
s~^, characteristic of the dSphs observed with MMFS. 
We further assume that the dSph has a Plummer sur- 
face brightness profile, such that members have projected 
radii drawn with probability p{R) oc (1 + E? /r^)~^ and 
sampled to a maximum radius of i? < bvp. We draw mag- 
nesium indices for the dSph members directly from the 
real MMFS measurements for the most likely {Pm > 0.9) 
Carina, Sculptor and Sextans members (we exclude For- 
nax, the most metal-rich of the observed dSphs, because 
of its systematically larger Mg values). 

We draw a fraction 1 — fmem of the simulated data 
points from a contaminant population that we assume 
to follow a uniform spatial distribution over the sampled 
region. We draw contaminant velocities from the Be- 
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TABLE 1 

Results from EM algorithm for MMFS data 



Galaxy 


N 


Nm 


( V^) mem 

(km s"-"-) 


^Vq ,raem 

(km s"-') 


(W>mem 

(A) 


(A) 


(A) 


(A) 


(J anna 
Fornax 
Sculptor 
Sextans 


1982 
2633 
1541 
947 


774 

2483 
1365 
441 


222.9 ±U.l 
55.2 ±0.1 
111.4 ±0.1 
224.3 ±0.1 


6.6 ± 1.2 
11.7±0.9 
9.2 ± 1.1 
7.9 ± 1.3 


0.40 ±0.01 
0.59 ± 0.01 
0.39 ± 0.01 
0.36 ± 0.01 


0.07 ±0.02 
0.12 ±0.02 
0.07 ±0.01 
0.05 ±0.02 


0.83 ±0.01 
0.51 ±0.01 
0.77 ±0.01 
0.72 ±0.01 


0.13 ±0.02 
0.26 ±0.06 
0.22 ±0.05 
0.21 ±0.03 




EllipUtal RftCljua (art5min> Hllptleai Badiu? (srcmin) 



Fig. 7. — Velocity distribution and membership fraction versus elliptical radius. The top sub-panels for each galaxy plot stellar velocity versus 
elliptical radius. Marker colors indicate Pm as in Figure^ The bottom sub-panels plot the membership fraction versus elliptical radius, estimated 
via monotonic regression. Dotted horizontal lines enclose the conventional 3cr velocity range. 



sangon model of the Carina foreground (Equation [71 see 
also the bottom-left panel of Figure[T]). We draw contam- 
inant EMg values from the real MMFS measurements of 
the least likely (Pm < 0.1) Carina, Sculptor and Sextans 
members. To each artificial data point we add a small 
scatter (in velocity and SMg) that corresponds to mea- 
surement errors drawn randomly from the real MMFS 
data. 

We consider the performance of the EM algorithm 
while varying the following three quantities: sample size 



N, member fraction fmem, and member mean veloc- 
ity {V)mem- The mean velocity controls the amount 
of overlap between member and contaminant distribu- 
tions. We consider three sample sizes to represent realis- 
tic dSph data sets. The smallest have 30 stars, typical of 
the sa mples now avai l able f o r the faintest nearb y dSphs 
(e.g., iSimon fc Gehal [20071 : iMartin et al.] [200l . The 
intermediate-size samples have 300 stars and are thus 
similar to many of the samples pu blished for brighter 
dSphs within the past five years (e.g.. lKlevna et a"Lll2002l : 
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Fig. 8. — Maps indicating probability of membership. As in Figure [Tl 
Black (red; magenta; green; cyan; blue) markers signify Pm ^ 0.01 {Pm > 
radii (the King outer radius of Sextans is beyond the window dimensions) 



Walker et al 



Mateo et al. 



[20061 : iMuiioz et al.l[2(M iKoch et al.l[2n07l : 
20081 ). The largest simulated samples have 



3000 stars, of the same magnitude as the MMFS samples. 
For each sample size, we consider three possible mem- 
ber fractions — fmem = 0-2,0.5,0.8 — that range from 
contaminant-dominated to member-dominated. For each 
case of sample size and member fraction, we then con- 
sider three possible values for the dSph mean velocity: 
{V)mem = 50, 100, and 200 km s^^. Since the artificial 
foreground peaks at a velocity of ^ 50 km s~^, these val- 
ues represent a range in degree of overlap between dSph 
and contaminant velocities. 

We apply the EM algorithm to each simulated data set 
and record the resulting estimates of the member mean 
velocity and velocity dispersion, as well as the member- 
ship probabilities determined for each star. Figures[91-IT71 
in Appendix [C] display the artificial data and demon- 
strate the effectiveness with which the EM algorithm dis- 
tinguishes members from contaminants in each of the 27 
simulations. Left-hand panels plot projected radius and 



marker color indicates probability, Pm , that the star is a dSph member. 
0.01; > 0.50; > 0.68; > 0.95; > 0.99). EUipses mark King core and "tidal" 
, adopted from IH95. 

EMg against velocity, using blue circles to represent sim- 
ulated dSph members and black squares to represent con- 
taminants. Middle panels show the same scatter plots, 
but with marker color now indicating, as in Figures [1] 
and [21 the membership probabilities returned by the EM 
algorithm. Right-hand panels identify any stars that are 
mis-classified by the EM algorithm — i.e., member stars 
for which the EM algorithm gives Pm < 0.5 and/or con- 
taminants for which Pm > 0.5. 

In 26 of the 27 simulations, the number of mis-classified 
stars amounts to less than 10% of the total sample, 
and the typical number of mis-classified stars is between 
1% — 3% of the total sample. In 23 of the 27 simulations, 
the EM algorithm returns accurate (i.e., the error bar 
includes the true value) estimates of the dSph mean ve- 
locity and velocity dispersion. Among the intermediate- 
sized and large samples {N > 300), the algorithm's lone 
failure occurs for the heavily contaminated {fmem — 0.2) 
N = 300 sample for which the dSph mean velocity of 
{V)mem = 50 km s^^ coiucidcs with the peak of the con- 
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taminant distribution. Under those circumstances the 
algorithm mis-classifies nearly two-thirds of all stars, and 
the measured velocity dispersion is effectively that of the 
foreground (see Figure fT2l bottom-left). 

The three other cases of inaccurate estimates all occur 
for simulations of small {N = 30) samples. As one might 
expect, the case with with the heaviest contamination 
ifmem = 0.2) and most velocity overlap {{V)mem — 50 
km s~^) yields an inflated estimate of the velocity dis- 
persion (see Figure [HI bottom left). The two remaining 
cases of failure at = 30 each merit some discussion. 
The simulation with heavy contamination (/mem = 0.2) 
and moderate overlap {{V)mem — 100 km s^^) yields an 
under-estimate of the velocity dispersion. This partic- 
ular failure results from a mis-classification of just two 
stars (among just six members). Inspection of Figure [TUl 
(bottom left) reveals that the member that is mistakenly 
classified as a contaminant occupies the extreme tails in 
the member distributions of velocity, position and Mg, 
while the contaminant mistakenly identified as a mem- 
ber lies near the center of the member distributions. A 
reasonable person examining the data shown in Figure 
[TOl (bottom left) would be likely to mis-classify these two 
stars in the same way as the algorithm. In the simulation 
with moderate contamination (/mem = 0.5) and minimal 
overlap {{V)mem = 200 km s~^) the algorithm correctly 
classifies all stars but, as happens with non-negligible 
probability for data sets of A^ < 15, the sample simply 
exhibits a smaller velocity dispersion than the population 
from which it is drawn (see Figure (TTJ top right). Thus 
these final two cases of failure do little to incriminate 
the EM algorithm, but rather remind us of the perils of 
working with small samples. We conclude that, barring 
a conspiracy among heavy contamination, close velocity 
overlap and small-number statistics, the EM algorithm 
provides an excellent tool for quantifying the properties 
of dSph kinematic samples. 

6.1. On the Faintest dSphs 

Mining of data from the Sloan Digital Sky Survey 
(SDSS) has recently uncovered a new population of ex- 
tremely faint dSphs in the Local Grou p, with absolute 
magn itudes as small as My ~ —2 (e.g.. iBelokurov et al.l 
|2007| ) . Initial kinematic studies of these systems suggest 
that they have velocity dispersions as small as ~ 4 km 
s""'^, systematically co lder than those of brighter dSphs 
(| Simon fc Gehal [20071 : [Martin et al.ll2007l[ ). In order to 
test the efficacy of the EM algorithm in the regime rel- 
evant to these systems, we consider nine additional sim- 
ulations that use samples of A^ = 30 stars drawn from 
member distributions with velocity dispersions of 4 km 
s^^. We again consider member fractions of fmem = 
0.2, 0.5, 0.8 and mean velocities (F)mem = 50, 100, 200 
km s^^). 

The results of these nine simulations (see Figures 
[IHH^ni of Appendix [Cj are similar to those for the small 
samples previously discussed. The EM algorithm re- 
turns accurate estimates in six of the nine simulations, 
the same rate of success as in the nine previous simu- 
lations with A^ = 30 and cvo^mem = 10 km s~^. Two 
of the three failures occur for the samples that have the 
heaviest contamination (fmem = 0.2) and most velocity 
overlap {{V) = 50, 100 km s^^). The third {fmem. = 0.5, 
{V)mem = 100 km s~^) barely fails, returning a veloc- 



ity dispersion of crvo,mem = 8.5 ± 4.0 km s~^ after mis- 
classifying just one star. With these simulations we find 
that the efficacy of the EM algorithm is insensitive to the 
member velocity dispersion, so long as that dispersion is 
sufficiently smaller — as is the case for all dSphs — than 
that of the contaminant population. Of course the cold- 
est systems require extreme caution if the velocity dis- 
persion is similar to the veloc ity errors, as is the cas e for 
some dSphs i n the s amples of ISimon fc Gehal (|2007f l and 
[Martin et al.l ([2007^). In such situations the algorithm's 
performance depends on external factors affecting the va- 
lidity of the velocity dispersion estimate itself, such as the 
degree to which the measurement errors are known (see 
Walker et al. in prep). 

6.2. Improvement via a velocity filter 

Because kinematic samples for the faintest dSphs are 
likely to suffer heavy contamination, one hopes to im- 
prove the reliability of the EM algorithm in that regime. 
Often it is possible to produce a reasonable guess of the 
mean velocity upon visual inspection of the V, R, and 
(if available) SMg scatter plots (e.g.. Figures [T8l--l20|. 
even for small and contaminated samples. In such cases, 
application of an initial velocity "filter" can restore the 
effectiveness of the algorithm. We have repeated all 36 
simulations after initializing all stars within 40 km s~^ 

of the dSph mean velocity to have PfP — 1, while all 

stars outside this range receive P^f = (recall that m 
the "unfiltered" algorithm, all stars receive initialization 
~ 0.5). We find that the filter corrects all noted 
cases of overestimated velocity dispersions. However, the 
filter does not correct the noted cases of underestimated 
dispersions for the A^ = 30 samples, as these failures 
continue to refiect the difficulties that arise from small- 
number statistics. Tables [2| and [3| list parenthetically the 
parameter estimates obtained after applying the velocity 
filter. 

7. COMPARISON WITH OTHER METHODS 

We now use the simulated data sets to compare pa- 
rameter and membership estimates obtained from the 
EM algorithm to those obtained using more conventional 
methods for contaminant removal. We consider two such 
methods. The first is the simple and widely used sigma- 
clipping algori thm. The s econd, advocated recently by 
[Klimcntowski e t al.l (I2007D and[Woitak & Lokas (2007), 
identifies and rejects stars that are not gravitationally 
bound to the dSph, according to a mass estimator based 
on the virial theorem. Like the EM algorithm, both 
methods require estimation of the mean velocity, and 
the sigma-clipping algorithm also requires estimation of 
the velocity dispersion. For consistency we continue to 
use Equations [T3| and [T4| to estimate these parameters. 
Here it is important to remember, however, one of the 
key differences between the EM and the alternative algo- 
rithms. While the EM algorithm deals with membership 
probabilities, the non-EM algorithms consider a given 
star as either a member or a nonmember. Therefore, in 
executing either of the non-EM algorithms, we initialize 
= 1 for all i, and stars that are rejected during 
the course of the algorithm receive Pm = during all 
subsequent iterations. 
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7.1. Sigma Clipping 

Sigma-clipping algorithms iteratively estimate the 
variance of some variable among member stars and then 
remove stars with outlying values. Historically, most 
dSph kinematic samples provide only velocity informa- 
tion (no line strength information). Hence a sigma- 
clipping routine typically involves estimating the mem- 
ber mean velocity and velocity dispersion, then reject- 
ing all stars with velocities that differ from the mean 
by more than some multiple of the velocity dispersion. 
The choice of threshold is arbitrary, with typical values 
between 2.5avo ^^'^ 4crvb- Here we consider a 3a clip- 
ping algorithm. In executing the (unfiltered) algorithm 

we initialize Pj^j^ = 1 for all i. Then each iteration of 
the algorithm consists of the following two steps: 1) es- 
timate the mean velocity and velocity dispersion via 50 
iterations of equations [T3l and [Ml then 2) assign Pm = 
to all stars with velocities more than 3ava from the mean. 
We terminate the algorithm after the iteration in which 
no stars are newly rejected. We record the estimates of 
the mean velocity and velocity dispersion obtained in the 
last iteration, and we estimate the number of members 
as X]t=i , which in this case equals the number of 
stars for which Pm = 1. 

7.2. Virial Theorem 



iKlimentowski et all (j2007D and IWoitak fc Lokad 
show that a second technique, which they apply 
after an initial sigma clip, can remove contaminants more 
reliably than sigma-clipping alon e. This method, origi- 
nally proposed in an anlaysis bv Iden Hartog fc KatgertI 
(|1996) of galaxy cluster kinematics, seeks explicitly to 
distinguish stars that are gravitationally bound to the 
dSph from those that are unbound. It incorporates 
a dynamical analysis that employs a mass estimator 
derived from the virial theorem. After sorting stars in 
sequence from smallest to largest projected radius and 
removing nonmcmbcrs (identified here as stars for which 

Pm = 0), the mass estimator is given by (jHeisler et al.l 



S7rNj:^^^{V,-{V)rnemf 



MvT{r) = 



J-1 



(16) 



where N is the number of stars with projected radius 
R < r and Dij is the projected distance between the 
i*'' and j*'' stars. Each rejection iteration is preceded by 
estimation of the mean velocity (and velocity dispersion, 
although this quantity is not used by the algorithm) us- 
ing Equations [13] and [TH In this "virial theorem" (VT) 
technique, a star at radius R is rejected (and assigned 
Pm = 0) if its velocity differs from {V)mem by more 
than yj2GMvT{r)/r evaluated at r = i?. Thus the VT 
algorithm goes beyond the sigma-clipping algorithm in 
that it considers the projected positions of the stars as 
well as their velocities. As in the EM algorithm, marginal 
outliers are less likely to be considered members if they 
are projected large distances from the dSph center, but 
here the explicit reason is that such stars are less likely 
to be gravitationally bound to the dSph. 

The VT algorithm requires initialization of (V^)mem 
and the values Pmi for all * (effectively a velocity cut); 



we consider two possible routes. Firs t, in what we denote 
as the "VT30-" algorithm, we follow [Klimentowski et all 
(2007") in using the 3tT clipping algorithm to provide the 
initialization. Second, in what we call the "VT^m" al- 
gorithm, we intialize {V)mem with the value obtained 
from the EM algorithm, and then assign Pm = to all 
stars for which the EM algorithm returned Pm < 0.5, 
and Pm = 1 to all stars for which the EM algorithm 
returned Pm > 0.5. 

Each iteration of the VT algorithm then consists of 
the following three steps: 1) estimate the mean veloc- 
ity and velocity dispersion via 50 iterations of equations 
[lS]and[T4l then 2) calculate Myrir) from Equation [TBI 
and finally 3) assign Pm = to all stars with veloci- 
ties more than ^j2GMvT{r)/r from the mean. As with 
the sigma clipping algorithm, we terminate the VT al- 
gorithm after the iteration in which no stars are newly 
rejected, recording final estimates of the mean velocity 
and velocity dispersion as well as the number of mem- 
bers, Pm, ■ 

7.3. Alternative Versions of the EM algorithm 

Finally, we consider alternative versions of the EM al- 
gorithm. The full version that we advocate and have 
described above uses velocity, magnesium and position 
information. But not all kinematic data sets have asso- 
ciated line strength information, and one may reasonably 
wonder about the effect of considering the stellar posi- 
tions. In order to examine the performance of the EM 
algorithm in the absence of line-strength data, and in 
order to examine the relevance of the positional informa- 
tion, we consider versions of the EM algorithm that do 
not use these bits of information. 

7.3.1. EMv: Velocity Only 

We denote as "EMy" a version of the EM algorithm 
that considers only the velocity data. Execution of the 
algorithm differs from that of the full EM algorithm in 
just two respects. In the expectation step we evalu- 
ate Pmem{Vi) using the univariate Gaussian probability 
given by Equation [6l rather than the bivariate probabil- 
ity Pmem(Vi,Wi) givcu by Equation [5l Second, we do 
not perform the monotonic regression estimate of p{a); 
instead, after the maximization step in the n*'' iteration 
we simply set p^"^(ai) equal to the global member frac- 
tion, N-^Hf^^Pl^^ , for all i. This effectively removes 
any influence of position on the membership probabili- 
ties and parameter estimates. Thus the EMy algorithm 
operates only on velocity data, similar in that respect 
to a conventional sigma clipping algorithm. The EMy 
algorithm differs from the sigma clip in that it assigns 
membership probabilities rather than enforcing an arbi- 
trary threshold, and it uses the expected distribution of 
foreground velocities to help evaluate these probabilities. 

7.3.2. EMvr: Velocity and Position Only 

We denote as "EMyi?" a version of the EM algorithm 
that considers the velocity and position data but ignores 
the magnesium index. This method is therefore applica- 
ble to most kinematic data sets. Its execution differs from 
that of the full EM algorithm only in the expectation 
step. There, we evaluate Pmem(^i) using the univariate 
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Gaussian probability given by Equation [6l rather than 
the bivariate probabihty Pmem{Vi,Wi) given by Equa- 
tion [5] 

7.3.3. EMvw: Velocity and Magnesium Only 

FinaUy, we denote as "FiMvw" a version of the EM al- 
gorithm that considers the velocity and magnesium data 
but ignores the stellar positions. By comparing its effec- 
tiveness to that of the full EM algorithm we can judge 
the the importance of including the position data. The 
execution of the EMyvK algorithm differs from that of 
the full EM algorithm only in that we do not perform 
the monotonic regression estimate of p(a). Instead, as in 
the EMy algorithm, after the maximization step in the 
nth iteration we set p^^^{ai) equal to the global member 
fraction, N-^E^^^Pj^^ , for all i, removing the influence 
of position on the membership probabilities and param- 
eter estimates. 

7.4. Comparison of Results Using Simulated Data 

Tables [2] (for the 27 simulations in which crvo.mem — 
10 km s~^) and [3] (for the nine simulations in which 
o'Vo,mem = 4 km s~^) list the results of applying 
each of the algorithms described above to the artifi- 
cial data sets. Underneath the row that identifies the 
data set, columns list the estimated number of mem- 
bers, Nmem = , followed by the number of mis- 
classified stars — i.e., the number of members for which 
Pm < 0.5 plus the number of nonmembers for which 
Pm > 0.5. Columns 5 and 6 list estimates of the dSph 
mean velocity and velocity dispersion, obtained from it- 
erating Equations 1 131 and 1 141 using the membership prob- 
abilities obtained from the algorithm. Parenthetical val- 
ues indicate estimates obtained after implementing the 
initial velocity filter described in Section [621 

We judge the performance of each algorithm accord- 
ing to whether or not it recovers (within the error bars) 
the mean velocity and velocity dispersion of the simu- 
lated member population. The final columns in Tables 
[2] and [3] indicate with either a "Y" or "N" whether the 
algorithm satisfies this criterion. As previously discussed 
(Section [H]) , the unfiltered EM algorithm succeeds in 23 
of 27 simulations with (Tv„,mem = 10 km s~^ (and in 6 
of 9 with avQ,mem = 4 km s~^), and performs well in all 
simulations in which the sample has either > 30 or 
fmem > 0.2. Usc of an initial velocity filter helps the al- 
gorithm's performance in cases of heavy contamination — 
the filtered EM algorithm succeeds in 32 of 36 simula- 
tions. 

7.4.1. EM versus sigma clipping and virial theorem 

In contrast, the unfiltered 3a algorithm succeeds in 
only four of 27 simulations with cryg^rnem = 10 km s-\ 
and in only two of nine with crvQ,mem = 4 km s~^. The 
unfiltered 3cr algorithm clearly has difficulty in the pres- 
ence of moderate to heavy contamination, failing in ev- 
ery simulation in which fmem < 0.8. Because the 3ct 
algorithm uses only velocity information, its reliability 
depends critically on the contrast between member and 
contaminant velocity distributions. In all cases of its fail- 
ure, the 3cr algorithm errs by failing to classify nonmem- 
bers as such, resulting in overestimates of the velocity 
dispersion. 



The 3(7 algorithm performs more reliably if we supply 
an initial velocity filter (see Section [6?2|) . The filtered 3cr 
algorithm succeeds in nine simulations for which it oth- 
erwise fails. However, we note that even without a filter, 
the EM algorithm generally performs well and distin- 
guishes members from contaminants more reliably than 
even the filtered 3cr algorithm. Clearly, one benefits by 
considering the contaminant distribution and using the 
additional information from positions and/or spectral in- 
dices when it is available. 

The performance of the YT^^ algorithm generally de- 
pends on the performance of the 3a algorithm that pro- 
vides its initialization. Counting all 72 filtered and unfil- 
tered runs, the VTs^ algorithm succeeds in all 21 cases 
in which the 3a algorithm succeeds, and fails in 43 of 
the 51 cases in which the 3a algorithm fails. In most 
cases the VT^^ algorithm manages to reject a handful 
of nonmembers that are unrecognized as such by the 3a 
algorithm, and so it yields slightly less inaccurate param- 
eter estimates. In eight cases the VT30. algorithm is able 
to compensate sufficiently to recover accurate paramter 
estimates where the 3a algorithm fails. Similarly, the 
performance of the VTem algorithm traces that of the 
EM algorithm, succeeding in all but one case in which 
the EM algorithm succeeds, and failing in all but one 
case in which the EM algorithm fails. 

It is important to note that the VT method has a 
different specific purpose than either of the methods 
chosen for its initialization. It seeks explicitly to iso- 
late a clean sample of gravitationally bound dSph mem- 
bers (Klimcntowski ct al. 2007; Wojtak & Lokas 2007). 
If such a sample is desired, as it is when applying equi- 
librium models, then the VT is appropriate to use after 
identifying likely foreground stars using the EM algo- 
rithm. The VT method is then complementary to the 
EM algorithm in this application. 

7.4.2. Which EM algorithm? 

The EM algorithm clearly outperforms the sigma clip- 
ping algorithm, and provides the appropriate initializa- 
tion for the VT algorithm when desired. But, assum- 
ing Mg data are available, which EM algorithm is best? 
Interestingly, the EMy^ algorithm, which excludes the 
magnesium information, enjoys the same rate of success 
as the full EM algorithm. The EMy^^ algorithm returns 
accurate parameter estimates in 23 of the 27 simulations 
with avo.mem = 10 km (and in six of the nine simula- 
tions with avo.mem = 4 km s~^). The EMv and EMy^y 
algorithms both succeed in 15 of 27 simulations with 
"■Vo.mem = 10 km (and six of the nine simulations 
with aVo,mem = 4 km s"^). 

These results prompt two conclusions. First, the supe- 
rior performance of the full EM and EMyji algorithms 
imply that it is important to consider the positions of 
stars when identifying contaminants. Given that the 
projected density of members declines as the density of 
contaminants remains approximately uniform, the a pri- 
ori probability of membership decreases with projected 
distance from the dSph center. Algorithms that incor- 
porate this trend outperform those that do not. Second, 
the ability of the EMvr algorithm to recover the correct 
parameters in some cases where the full EM algorithm 
fails suggests that there is room to improve our model 
of the Mg index distribution. Where the full EM algo- 
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rithm fails, it fails because it mis-classifies contaminants 
with small Mg-index values as members. These con- 
taminants correspond to metal-poor giants in the Milky 
Way halo, whereas our Gaussian model of the contami- 
nant Mg distribution allows for only a single component 
that is dominated by metal-rich dwarfs in the Milky Way 
disk. The full EM algorithm may therefore benefit from 
a contaminant model in which the Mg distribution is 
a double-Gaussian, allowing for a distinct halo compo- 
nent. We have not adopted such a model here because 
the single-component model performs well in all simu- 
lated data sets of comparable size to our MMFS samples. 
A two-component foreground model may become neces- 
sary, however, to optimize performance of the EM algo- 
rithm with small, heavily contaminated samples typical 
of those obtained for ul tra-faint dSphs ()Simon fc Gihil 
[2007t iMartin et al.|[2007h . 

8. DISCUSSION AND SUMMARY 

We have introduced the EM method as an effective tool 
with which to attack the problem of sample contamina- 
tion. In our implementation, the EM algorithm operates 
on all the available velocity, magnesium-index and po- 
sitional information. This marks a significant improve- 
ment over contaminant-removal methods that rely on ve- 
locity alone, particularly the conventional sigma-clipping 
method. Further, with appropriate modification of the 
likelihood function (Equation [T]), the EM algorithm can 
easily be tailored to incorporate other sorts of data as 
diagnostics of membership (e.g., photometric data, as in 
[Gilbert et a"n[20061 Apphed to our MMFS data, the EM 
algorithm evaluates the probability that a given star is a 
dSph member, while simultaneously estimating the mean 



velocity, velocity dispersion, and fraction of members as 
a function of position. 

Withi n the context of app l ying e quili brium dynamical 
model s. iKlimentowski et ail ()2007f ) and lWoitak fc Lokai 
discuss the importance of removing a second type 
of contaminant — former dSph members that have been 
tidally stripped and are thus no longer gravitationally 
bound to the dSph. These authors prescribe an itera- 
tive procedure that uses the virial theorem to provide a 
nonparameteric estimate of the dSph mass profile, and 
then removes stars with velocities greater than the local 
circular velocity (Section [721) ■ For the following reasons, 
this method should not be viewed as an "alternative" to 
the EM algorithm described above. First, while it uses 
velocity and position in applying the virial theorem, it 
does not use SMg, which we have shown to be a pow- 
erful, independent indicator of dSph membership. Sec- 
ond, and more fundamentally, it has a different and more 
specific purpose — the identification of a clean sample of 
gravitationally bound dSph members. If such a sample is 
desired, as it is when applyin g equilibrium model s , then 
the virial theorem me t hod o f IKlimentowski et al.l ()2007l ) 
and IWojtak fc LokasI ()2007t ) can still be applied appro- 
priately after identifying likely foreground stars using the 
EM algorithm. 
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uate School at the University of Michigan for generous 
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servatory. MM acknowledges support from NSF grants 
AST-0206081 0507453, and 0808043. EG acknowledges 
support from NSF Grants AST-0205790, 0505711, and 
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TABLE 2 

Results of Algorithms'' with Simulated Data and av^ 



10 KM s 



N 


^mem 


^wrong 








Success? 








(km s-i) 




(km s-i) 





Simulation 

EM 

EMv 

EM\/i4/ 
So- 

VTbm 

Simulation 

EM 

EMv 
EMvH 
EMvw 
So- 

VTgM 

Simulation 

EM 

EMv 
EMvfl 
EMvw 
So 

Simulation 

EM 

EMv 

EMvi; 

EMviv 

So 

VTs,, 

Simulation 

EM 

EMv 

EMvfl 

EMviv 

So 

Simulation 

EM 

EMv 

EMv« 

EMvw 

So- 

Simulation 

EM 

EMv 

EMvfl 

EMvw 

So- 

VTbm 

Simulation 

EM 

EMv 

EMvfl 

EMvw 

So 

Simulation 

EM 

EMv 

EMvfl 

EMvw 

So- 

VTgM 



N = 30 



N = 30 



N = 30 



N = 30 



N = 30 



N = 30 



N = 30 



N = 30 



N = 30 



Nmem — 

23(23) 
23(23) 
23(23) 
23(23) 
29(27) 
26(25) 
22(22) 



24 



15 



16(16) 
12(12) 
15(15) 
16(16) 
30(24) 
28(21) 
13(13) 

Nmem — 6 

8(5) 
13(6) 
14(5) 
12(4) 
29(16) 
28(9) 

6(4) 



N 



mem — 

23(23 
23(23 
24(24 
23(23 
30(27" 
27(25 
22(22 



24 



15 



15(15 
15(15 
15(15 
14(14 
30(21 
30(15 
14(l4 



5(5) 

23(8) 
11(11) 

5(5) 
30(17' 
30(11 

5(5) 



N 



N 



mem — 


23 


23) 


24 


24j 


24 


24) 


23 


23) 


30 


24) 


28 


22) 


23 


23) 


mem — 


14 




14 


\% 


14 


14) 


14 


14) 


30 


15) 


30 


is) 



24 



15 



13(13) 



5(5) 
13(5) 
14(5) 

9(5), 
30(6 
30(3 









2 


2 


1 


1 








5 


3 


2 


3 


2 


2 



2(2 

5(5 

5(5 

2(2 

15(9 

13 8' 

2(2) 



2(0) 
20(5) 
10(0) 

5 1) 
23(10) 
22(3) 

0(2) 









1 


1 


1 


1 








6 


3 


00 


1 


2 


2 



1(1) 

5(5) 
1(1) 

2(2) 
15(6 
15(2 
3(3) 



2(2) 
24(10) 
8(8) 

24(11) 
24(5) 
3(3) 









1 


1 


1 


1 








6 


2 


4 


2 


1 


1 



0(0 
0(0 
0(0 
0(0 
15(0 
15(2 
2(2) 



0(0 
4(0 
7(0 
4(0 
24(0 
24 3' 



(V)mem - 

48.8 ± 0.6( 
47.9±0.5( 
48.9±0.6( 
48.7±0.6( 
51.0±0.5( 
48.5±0.5( 
48.6 ± 0.6( 

(V)mem - 

46.1±0.5( 
48.2±0.5( 
46.3±0.4( 
46.2±0.5( 
48.8±0.5( 
48.8 ± 0.5( 
45.1 ±0.4( 



= 50kms~^ 
'48.8 ± 0.6) 
"47.8 ± 0.5) 
'48.9 ± 0.6) 
'48.7 ±0.6) 
'47.0 ± 0.5) 
'46.1 it 0.5) 
|48.6 ± 0.6) 

= 50kms~^ 
'46.1 ±0.5) 
'48.2 ± 0.5) 
"46.3 ± 0.4) 
'46.2 ± 0.5) 
'47.6 ± 0.4) 
'46.8 ± 0.4) 
'45.1 ±0.4) 



(V)mem = 50kms-^ 
34.5 ± 1.3(49.2 ± 1.0) 
20.2 ±0.9(43.5 ±0.9) 
28.5 ± 0.9(49.4 ± 1.1) 
32.0 ± 1.0(48.5 ± 0.8) 
45.0 ±0.7(40.1 ±0.8) 
42.2 ±0.8(43.8 ±1.0) 
49.2 ±1.0(51.2 ±1.4) 

{V)„,em = lOOkms-l 
102.6 ±0.5(102.6 ±0.5 
101.8 ±0.5(101.8 ±0.5 
102.3 ±0.5(102.3 ±0.5 
102.6 ±0.5(102.6 ±0.5 

93.4 ±0.5(99.8 ±0.5) 
99.8 ±0.5(101.2 ±0.5) 
102.8 ±0.5(102.8 ±0.5) 



(V) mem 
102.9 ±0 

99.1 ±0 
102.7 ±0 
101.5 ±0 

75.3 ±0 

75.3 ± 
102.7 ± 

(V) mem 
89.6 ±0 

85.2 ±0 
104.2 ± 

90.6 ± 

81.4 ± a 
81.4 ± a 
87.8 ± 



= lOOkms"^ 
6(102.9 ±0.6) 
6(99.1 ± 0.6) 
6(102.7 ±0.6 
6(101.5 ±0.6 
5(97.7 ±0.5) 
5(103.0 ±0.6) 
6(102.7 ±0.6) 

= lOOkms"^ 
9(89.6 ± 0.9) 
6(105.8 ±0.8) 
6(104.3 ±0.7) 
9(90.6 ± 0.9) 
6(102.1 ±0.7) 
6(104.6 ±0.8) 
9(87.8 ± 0.9) 



(V) mem — 
199.7 ±0.6 
199.6 ±0.6 

199.6 ±0.6 

199.7 ±0.6 

174.8 ± 0.6 

184.9 ± 0.7 
199.3 ±0.6 



200kms"^ 
199.7 ±0.6 
199.6 ±0.6 

199.6 ±0.6 

199.7 ±0.6 

200.2 ± 0.6 
199.9 ± 0.6 

199.3 ±0.6 



(V)„,em = 200kms~^ 
202.7 ±0.4(202.7 ±0.4 

202.7 ±0.4(202.7 ±0.4 
202.5 ±0.4(202.5 ±0.4) 

202.8 ± 0.4(202.8 ± 0.4) 

123.9 ± 0.5(202.7 ± 0.4) 
123.9 ±0.5(201.9 ±0.4) 
201.9 ±0.4(201.9 ±0.4) 



(V) mem 

198.2 ±0 
125.8 ±0 

128.3 ±0 
154.3 ± 
89.7 ±0 
89.7 ±0. 
201.7 ± 1 



= 200kms~^ 
7(198.2 ±0.7" 
5(197.9 ±0.7' 
6(198.5 ± 0.7' 
7(197.6 ± 0.7' 
5(198.3 ±0.7 
)(201.7± 1.1' 
1(201.7 ± 1.1) 




= lOkms"^ 
12.2 ±3.2) 

13.7 ±3.0) 

12.8 ± S.S) 

12.2 ± 3.2) 

15.3 ± 3.2) 
15.1 ±3.3) 
12.6 ±3.3) 



7.9±2.7i 

5.3 ± 1.9i 
7.4±2.6i 
7.9±2.7i 
42.6±6.5i 
28.5±6.4i 
8.0±2.6i 



lOkms" 



24.2 ±8.1(5.1 ±3.1) 
29.7 ±8.1(7.9 ±3.2) 

22.3 ±6.5(5.2 ±3.1) 
27.0 ±7.4(4.7 ±2.5) 
48.7 ±8.4(13.5 ±4.3) 
47.3 ±8.2(8.9 ±3.9) 
5.1 ±3.0(5.8 ±3.7) 

OVo mem ~ lOkmS""*^ 

7.8'± 2.7(7.8 ±2.7 

6.7 ±2.3(6.7 ±2.3' 
7.5 ±2.5(7.5 ±2.5 

7.8 ±2.7(7.8 ±2.7' 
22.3 ±4.0(11.7 ±2.9 
11.7 ±2.9(10.0 ±2.9 

8.0 ±2.7(8.0 ±2.7) 



o-v, 
10.5 ±3.3 

10.5 ±3.2 
10.8 ±3.3 

9.6 ±2.9 

38.6 ±6.2 
38.6 ±6.2 
11.1 ±3.5 



= lOkms"^ 
10.5 ± 3.3 
10.5 ± 3.2 
10.8 ± 3.3) 
9.6 ±2.9) 

14.0 ± 3.9) 
10.8 ± 3.4) 

11.1 ± 3.5) 



O-Vo mem ~ lOkmS 

5.S'± 2^7(5.3 ± 2.7) 

43.1 ±7.2(18.3 ±4.6) 
18.3 ±4.6(18.0 ±4.6) 

4.7 ±2.6(4.7 ±2.6) 
44.9 ± 7.8(20.6 ± 5.2 
44.9 ±7.8(18.3 ±5.0' 

4.2 ± 2.6(4.2 ± 2.6) 

O-Vo mem ~ 10kmS~- 

8.6'±T2(8.6 ± 3.2 

8.4 ±3.1(8.4 ±3.1' 

8.5 ±3.1(8.5 ±3.1' 

8.6 ±3.2(8.6 ±3.2' 
57.3 ±6.4(8.1 ±2.9 

44.2 ±5.1(8.1 ±2.8 
8.5 ±3.1(8.5 ±3.1) 

O-Vo mem ~ lOkmS"-"- 

4.5'± L8(4.5 ± 1.8 

4.4 ± 1.8(4.4 ± 1.8' 
4.1 ± 1.7(4.4 ± 1.7 

4.5 ±1.8(4.5 ±1.8' 
86.8 ± 9.6(4.5 ± 1.8 
86.8 ± 9.6(4.0 ± 1.4 
4.0 ±1.4(4.0 ±1.4) 

O-Vo mem ~ lOkmS^"' 

9.6'±r2(9.6 ±4.2) 
88.2 ±9.4(9.1 ±4.0) 
90.8 ±11.1(9.6 ±4.2) 
60.0 ±8.3(9.1 ±3.9' 
81.7 ±8.7(9.7 ±4.2 
81.7 ±8.7(12.0 ±5.5) 
12.0 ±5.5(12.0 ±5.5) 



Y(Y) 
N(N) 

y(y) 

Y(Y) 
N(N) 
N(N) 
Y(Y) 



= lOkms" 


1 






(7.9 ±2.7 




Y 


Y) 


5.3 ± 1.9 




N 


N) 


(7.4 ±2.6 




Y 


y) 


7.9 ± 2.7 




Y 


Y) 


15.9 ± 3.4) 


N 


N) 


13.2 ± 3.1) 


N 


N) 


(8.0 ± 2.6) 


Y 


y) 



N(N) 


N 


Y) 


N 


N) 


N 


N) 


N 


Y) 


N 


Y) 


N 


n) 



Y 




N 


I 


Y 


Y 


Y 


Y 


N 


Y 


N 


Y 


Y 


Y 



Y(Y 
Y(Y 
Y(Y 
Y(Y 
N(N 
N(Y 
Y(Y 



N 


N 


N 


N 


N 


N 


N 


N 


N 


N 


N 


N 


N 


N 



Y 


Y 


Y 


Y 


Y 


Y 


Y 


Y 


N 


Y 


N 


Y 


Y 


Y 



N 


N 


N 


N 


N 


N 


N 


N 


N 


N 


N 


N 


N 


N 



Y 




N 


?! 


N 


Y) 


N 


Y) 


N 


Y) 


N 


Y) 


Y(Y) 



EM Algorithm 
TABLE 2 — Continued 



N 


^mem 


^wrong 






Success? 








(km s-l) 


(km s ) 





Simulation 

EM 

EMv 

EMvfl 

^Mvw 

3a 

Simulation 

EM 

EMy 

EMvR 
EMvw 
3a- 

VTs,, 
VTbm 

Simulation 

EM 

EMy 

EMvR 

EMvw 
3a 



N = 300 



N 



mem 


= 240 






(V)mem - 


= 50kms"^ 


242(241 




6(6) 


48.7 ±0.2(48.7 ±0.2) 


235 


234 




25 


25) 


48.7 ±0.2 


48.7 ±0.2) 


238 


238 




15 


15) 


48.6 ±0.2 


48.6 ± 0.2) 


243 


241 




9 


8) 


48.9 ±0.2 


48.7 ±0.2) 


271 


271, 




31 


31) 


48.4 ±0.2 


48.4 ± 0.2) 


255 


255* 




23 


23) 


48.5 ±0.2 


48.5 ± 0.2) 


223(223 




21 


21) 


48.7 ±0.2(48.7 ±0.2) 



N = 300 



N = 300 



Nmem = 150 

153(153) 
149(148) 
147(147) 
156(152) 
273(238) 
223(195) 
130(1,30) 

Nrnem = 60 

128(69) 
72(72) 
74(74) 
129(67) 

293(188) 
280(128) 
46(51) 



EM 

EMy 

EMyfl 

EMy^y 

3a 

Simulation 

EM 

EMy 
EMyH 

EMvw 
3a 

VTs, 

Simulation 

EM 

EMy 

EMvR 
EMvw 
3a 

Simulation 

EM 

EMy 
EMyfl 

EMvw 
3a 

VTs,, 



67(67) 
79(70) 

78(78) 
120(65) 
294(169) 
280(105) 
54(54) 



N = 300 



N: 



N = 300 



N 



N = 300 



N: 



mem — 

60(60 
62(62 
60(60 
85(61 
298(73 
296(54' 
51(51) 



60 



10(10 
51(51 
37(37 
17(15 
123(88) 
79(53) 
26(26) 



197(11) 
66(66) 
21(21) 
199(17) 
233(128) 
220(68) 
106(13) 



Simulation 


N = 


300 


Nmem 


= 240 






EM 






241(241 




8(8) 


EMy 






244 


244 




23 




EMvR 






241 


241 




17 


Si 


EMy^y 






240 


240 




11 


11) 


3a 






261 


261 




23 


23) 


VT3a 






237 


237 




21 


21) 








219 


219 




21 


21) 


Simulation 


N = 


300 


-Nmem 


= 150 






EM 






151 


151 




9(9) 


EMy 






156 


155 




37(36) 


EMvR 






147 


147 




22(22) 


EMyvi^ 






155 


154 




16 


16) 


3a 






296 


216 




146(66) 








282 


173 




132(41) 








123 


123 




33(33) 


Simulation 


N = 


300 


-Nixiem — 


60 







9(9) 
56(58) 

18(19) 
60(18) 
234(109) 
220(45) 
12(12) 



mem = 240 






242 


242 


3 


3 




242 


242 


5 


5 




241 


241 


4 


4 




242 


242 


3 


3 




243 


243 


5 






225 


225 


23 


23) 


225 


225 


21 


21) 


mem — 150 






150 


150 


3 






149 


149 


5 


t 




151 


151 


3 






150 


150 


2 






300 


155 


150(5) 


299 


126 


149(24) 


124 


124 


26(26) 



9 

23(2) 
238(13) 
236(6) 

9(9) 



(V)mem = 50kms-^ 

50.3 ±0.2(50.3 ±0.2 

50.0 ±0.2(50.0 ±0.2 

50.1 ±0.2(50.1 ±0.2) 
50.1 ±0.2(50.4 ±0.2) 
52.9 ±0.2(50.8 ±0.2) 

51.4 ±0.2(50.4 ±0.2) 

50.1 ±0.3(50.1 ±0.3) 

(V>i„em = 50kms-^ 
70.9 ±0.2(50.5 ±0.3) 
48.4 ±0.2(48.4 ±0.2) 
50.4 ±0.3(50.4 ±0.3) 
70.9 ± 0.2(49.7 ± 0.3) 

63.3 ±0.2(49.2 ±0.2) 

63.4 ±0.2(50.7 ±0.2) 

67.2 ±0.3(52.1 ±0.4) 

(V)mem = lOOkmS-1 

99.6 ±0.2(99.6 ±0.2) 

99.7 ±0.2(99.7 ±0.2) 
99.7 ±0.2(99.7 ±0.2) 

99.5 ± 0.2(99.5 ± 0.2) 

99.3 ±0.2(99.3 ±0.2) 

99.4 ±0.2(99.4 ±0.2) 
99.2 ±0.2(99.2 ±0.2) 

lOOkms"^ 
99.0 ±0.2 

98.5 ± 0.2 
99.0 ± 0.2 
98.7 ±0.2 

96.6 ± 0.2 
97.6 ± 0.2 
99.6 ±0.2 



(V)mem — 

99.0 ±0.2 
98.4 ±0.2 

99.0 ±0.2 
98.6 ±0.2 
82.9 ±0.2 

84.1 ±0.2 

99.6 ±0.2 

(V)m 

98.0 ±0.3 

94.7 ±0.3 

97.2 ±0.3 
87.7 ±0.3 
76.7 ±0.2 
76.2 ±0.2 

99.1 ±0.4 

(V)mem — 

198.7 ±0.2i 
198.5 ±0.2 
198.7 ±0.2 
198.7 ±0.2 

198.5 ±0.2 
198.7 ±0.2 
198.9 ±0.2 

(V) mem — 

198.7 ±0.2 

198.6 ±0.2 
198.4 ± 0.2 

198.8 ± 0.2 
134.6 ± 0.2 
133.8 ± 0.2 
198.8 ±0.2 



lOOkms-^ 
98.0 ± 0.3) 

95.8 ±0.3) 
97.2 ± 0.3) 

97.9 ± 0.3) 
93.2 ± 0.2) 

94.0 ± 0.3 

99.1 ±0.4; 

200kms~^ 
198.7 ± 0.2) 
198.5 ± 0.2) 

198.7 ±0.2) 
198.7 ±0.2) 

198.5 ±0.2) 
198.7 ±0.2) 
198.9 ±0.2) 

200kms~^ 

198.7 ±0.2 

198.6 ±0.2 

198.4 ±0.2 

198.8 ±0.2 

198.5 ±0.2 
198.8 ±0.2 
198.8 ±0.2 



{V)mem = 200kms-^ 

199.5 ±0.3(199.5 ±0.3 
198.1 ±0.3(198.1 ±0.3 
199.7 ±0.3(199.7 ±0.3 
168.4 ±0.3(199.0 ±0.3 
95.9 ±0.2(196.4 ±0.3 
96.4 ±0.2(199.3 ±0.3 
199.3 ±0.3(199.3 ±0.3) 



(TVo mem ~ lOkmS 

10.8'±T9(10.7± 1.9) 
10.7 ±1.8(10.6 ±1.7) 

10.7 ±1.8(10.7 ±1.8) 
11.2 ± 1.9(10.8 ± 1.9) 

12.6 ±2.1(12.6 ±2.1) 
11.9 ±2.0(11.9 ±2.0) 

10.8 ±1.9(10.8 ±1.9) 

(TVo j„em ~ lOkmS"""" 

10.9 ±To(10.9 ± 2.0) 
9.7 ±1.7(9.7 ±1.7) 
9.2 ±1.8(9.2 ±1.8) 

12.2 ±2.1(11.4 ±2.0) 

24.7 ±2.9(16.1 ±2.2) 
16.9 ±2.3(13.5 ±2.1) 

10.1 ±2.0(10.1 ±2.0) 

(TVq n,e„ = lOkms"""" 
60.6'± 4.9(8.8 ± 2.3) 
8.2 ±1.6(8.2 ±1.6' 
8.6 ±2.0(8.6 ±2.0 
60.5 ±4.9(8.8 ±2.2) 

48.2 ±4.0(17.4 ±2.6" 

45.0 ±3.9(16.7 ±2.7 
48.9 ±6.1(8.9 ±2.6) 

asjg = lOkms"""" 
10.2'± L8(10.2 ± 1.8) 

10.1 ± 1.7(10.1 ± 1.7) 
9.9 ±1.7(9.9 ±1.7) 

10.2 ±1.8(10.2 ±1.8 
10.7 ±1.8(10.7 ±1.8 

9.9 ±1.7(9.9 ±1.7 
9.6 ±1.7(9.6 ±1.7' 



10.4± 1.9i 
11.9 ±1.9 
10.4 ± 1.8 
11.3 ±2.0 
33.3 ±3.4 
29.8 ±3.2 
10.1 ±2.1 



= lOkms" 
10.4 ± 1.9 
11.8 ± 1.9' 

10.4 ± 1.8 
11.3 ± 2.0 

15.5 ± 2.4' 
13.8 ± 2.3 

10.1 ± 2.1; 

= 10kms~- 



(TV, 

9.5 ± 1.8(9.5 ± 1.8 

9.3 ±1.8(9.3 ±1.8 
9.5 ± 1.8(9.5 ± 1.8 

9.4 ±1.8(9.4 ±1.8' 
82.3 ± 5.3(9.6 ± 1.9 
81.2 ± 5.4(9.7 ± 2.0 

9.8 ±2.0(9.8 ±2.0) 

(^Cmem = lOkmS-^ 

9.2 ±2.1(9.2 ±2.1) 
9.7 ±2.3(9.7 ±2.3) 

9.3 ±2.1(9.3 ±2.1) 
63.2 ± 7.0(9.5 ± 2.2) 
70.9 ±5.1(12.1 ±2.5) 
70.8 ±5.1(9.4 ±2.2) 

9.7 ±2.3(9.7 ±2.3) 



Y(Y) 


Y 


Y) 


Y 


Y) 


Y 


Y) 


N 


Y) 


Y 


Y) 


Y 


Y) 



Y(Y) 
N(Y) 
N(N) 
N(N) 
Y(Y) 

N(Y) 
N(N) 
Y(Y) 
N(Y) 
N(N) 
N(N) 
N(Y) 



Y(Y) 


Y 


Y) 


Y 


A 


Y 


Y) 


Y 


Y) 


Y 


Y) 


Y 


Y) 



) Y 


Y) 


) Y 


Y) 


Y 


Y) 


Y 


Y 


N 


N) 


N 




) Y 


?! 



10.7 ±2.5(10.7 ±2.5) Y(Y) 
12.9 ± 2.2(9.9 ± 1.9) N(Y) 

11.8 ±2.4(11.7 ±2.4) Y(Y) 
39.6 ±4.2(11.5 ±2.5) N(Y) 

48.9 ±4.2(18.6 ±2.8) N(N) 



10.6 ±2. 



(16.0 ± 2.9 


N 


N 


(10.6 ±2.8 


Y 


Y 


= lOkms"^ 




(10.2 ± 1.7 


Y 


Y 


(10.4 ±1.7 


Y 


Y 


(10.2 ± 1.7 


Y 


Y 


(10.2 ± 1.7 


Y 


Y 


(10.2 ± 1.7 


Y 


Y 


(10.1 ± 1.7 


Y 


Y 


(10.2 ± 1.7 


Y 


Y 


= lOkms"^ 





Y 


Y 


Y 


Y 


Y 


Y 


Y 


Y 


N 


Y 


N 


Y 


Y 


Y 



Y 


Y 


Y 


Y 


Y 


Y 


N 


Y 


N 


Y 


N 


Y 


Y 


Y 
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TABLE 2 



Continued 



N 




^wrong 








Success? 








(km s-l) 




(km s ) 





= lOkms"^ 




(10.1 ± 1.0 


Y 


Y 


(10.2 ± 0.9 


Y 


Y 


(10.1 ± 1.0 


Y 


Y 


(10.2 ± 1.0 


Y 


Y 


(11.2 ± 1.0 


N 


Y 


(10.4 ± 1.0 


Y 


Y 


(10.0 ± 1.0 


Y 


Y 


= lOkms"^ 




(10.6 ± 1.1 


Y 


Y 


(10.4 ± 1.0 


Y 


Y 


(10.2 ± 1.0 


Y 


Y 


(10.9 ± 1.1 


Y 


Y 


(15.0 ± 1.3 


N 


N 


(13.1 ± 1.2 


N 


N 


(10.3 ± 1.1 


Y 


Y 


= lOkms"^ 





Simulation N : 

EM 
EMv 

EMyfl 

^Mvw 

3(T 

VT3,. 

Simulation N : 

EM 

EMy 

FMvR 
'EUvw 

3tT 

VTs,, 

Simulation N : 

EM 

EMy 
EMyfl 

Simulation N : 

EM 

EMy 
EMyji 

3(T 

VTsa 

VTbm 

Simulation N : 

EM 

EMy 

EMyfl 

EMyvi^ 

30- 

VTsa 

VT_g;M 



3000 



3000 



3000 



3000 



3000 



Nmem = 2400 

2417(2416) 
2432(2430) 
2421(2421) 
2416(2413) 
2696(2696) 
2492(2492) 
2260(2260) 

Nmem = 1500 

1540(1538) 
1569(1565) 
1545(1544) 
1566(1555) 
2783(2340) 
2341(1927) 
1305(1305) 

Nmem = 600 

676(674) 
662(650) 
653(650) 
718(710) 
2894(1936) 
2719(1337) 
446(446) 

Nmem = 2400 

2403(2402) 
2410(2408) 
2406(2405) 
2399(2398) 
2599(2599) 
2399(2399) 
2226(2226) 

Nmem = 1500 

1516(1516) 
1522(1521) 
1528(1527) 
1520(1518) 
2917(2104) 
2630(1689) 
1283(1283) 



Simulation N : 

EM 

EMy 
EMyfl 

'EMvw 
3a 

VTsa 

VT^M 

Simulation N : 

EM 

EMy 
EMyH 

EMvw 
3a 

VT3a 



3000 



3000 



Nmem = 2400 

2402(2402) 
2393(2392) 
2393(2393) 
2402(2401) 
2401(2401) 
2231(2231) 
2237(2237) 

Nmem = 1500 

1501(1501) 
1502(1502) 
1500(1500) 
1503(1503) 
2994(1545) 
2978(1289) 
1280(1280) 



71(71) 
277(277) 
198(199) 
109(107) 
300(300) 
236(236) 
196(196) 



121(121) 

583(584) 
295(296) 
192(191) 
1283(840) 
849(499) 
253(253) 



134(134) 
581(587) 
241(241) 
210(209) 
2294(1336) 
2123(777) 
204(204) 



47(47) 
190(190) 
143(144) 

71(72) 
207(207) 
187(187) 
196(196) 



86(86) 
394(394) 
230(230) 
125(127) 
1417(604 
1140(303 
241(241) 



Simulation N = 3000 


Nmem = 600 




EM 


641(640) 


99(98) 


EMy 


649(620) 


447(439) 


EMvR 


601(599) 


198(198) 


EMyw 


933(655) 


362(154) 


3a 


2910(1570) 


2310(970) 


VTsa 


2804(955) 


2206(421) 




427(427) 


205(205) 



12 


12) 


24 


24) 


19 


19) 


10 


10) 


31 


31) 


175 


175 


169 


169 



15(15) 
57(57 
45(45) 
22(22) 
1494(55) 
1478(221) 
224(224) 



(V), 
50.2 ±0.11 
50.1 ±0.11 

50.1 ±0.11 

50.2 ±0.11 
50.1 ±0.11 

50.1 ±0.11 

50.2 ±0.11 

(V)mem - 

49.6 ±0.11 
49.9 ±0.11 

49.6 ±0.11 

49.8 ±0.11 

54.3 ±0.11 

51.7 ±0.11 

49.7 ±0.11 

(V)mem - 

50.1 ±0.11 
51.0 ±0.11 

50.6 ±0.11 

49.9 ±0.11 
62.5 ±0.11 

60.8 ±0.11 

50.5 ±0.11 

mem — 

100.2 ±0.1l 
100.2 ±0.1l 
100.2 ±0.11 
100.2 ±0.11 

99.8 ±0.11 

99.9 ±0.11 
100.2 ±0.11 

(V) mem — 

99.7 ±0.11 

99.7 ±0.11 

99.6 ±0.11 

99.6 ±0.11 

83.8 ±0.11 
87.0 ±0.11 
99.8 ±0.11 

(V) mem — 
99.5 ±0.1 

98.7 ±0.1 
100.1 ±0.1 

90.7 ±0.1 

73.4 ±0.1 

73.5 ±0.1 
100.1 ±0.1 

mem — 

200.0 ±0.1 

200.1 ±0.1 
200.0 ±0.1 
200.0 ±0.1 
200.0 ±0.1 
200.0 ±0.1 

200.0 ±0.1 

mem — 

199.9 ±0.1 

200.1 ±0.1 
200.0 ±0.1 
199.9 ±0.1 
136.6 ±0.1 
136.8 ±0.1 
200.0 ±0.1 



= 50kms~^ 

'50.2 ± 0.1) 
'50.1 ± O.l) 
'50.1 ± O.l) 
'50.2 ± 0.1) 
'50.1 ± O.l) 
'50.1 ± O.l) 
|50.2±0.l) 

= 50kms~^ 

'49.6 ±0.1) 
'49.9 ± O.l) 
'49.6 ± O.l) 
'49.7 ± O.l) 
'50.1 ± 0.1) 
'50.0 ± O.l) 
|49.7±0.l) 

= 50kms~^ 

'50.1 ±0.1) 
'51.0 ± O.l) 
'50.6 ± O.l) 
'50.0 ± O.l) 
'50.6 ±0.1) 
'50.5 ±0.1) 
|50.5±0.1) 

: lOOkms^i 
100.2 ±0.1) 
'100.2 ±0.1) 
'100.2 ±0.1) 
'100.2 ± O.l) 
'99.8 ±0.1) 
'99.9 ±0.1) 
|100.2±0.1) 

: lOOkms^^ 
'99.7 ±0.1) 
'99.7 ± O.l) 
'99.6 ± O.l) 
'99.7 ± O.l) 
'98.0 ±0.1) 
'98.6 ±0.1) 
|99.8±0.l) 

lOOkms"^ 
'99.5 ± 0.1) 
'99.0 ± O.l) 

100.2 ±0.1) 
'99.0 ±0.1) 
'95.3 ±0.1) 
'97.1 ±0.1) 
'l00.1±0.1) 

200kms"^ 
'200.0 ±0.1) 
'200.1 ±0.1) 
'200.0 ±0.1) 
'200.0 ± 0.1) 
'200.0 ± 0.1) 
'200.0 ± 0.1) 
'200.0 ± 0.1) 

200kms"^ 
199.9 ±0.1) 
'200.1 ±0.1) 
'200.0 ± 0.1) 
'199.9 ± 0.1) 
199.8 ± 0.1) 
'200.0 ± 0.1) 
'200.0 ± 0.1) 



o-Vn 



10.6 ±1.1 
10.5 ± 1.0 

10.2 ± 1.0 
11.1 ± 1.1 

27.3 ±1.7 

18.7 ± 1.4 
10.3 ± 1.1 



10.3 ±1.3(10.2 ±1.3) 
10.1 ± 1.0(9.8 ± 1.0) 
9.9 ±1.1(9.8 ±1.1) 
11.5 ± 1.3(11.1 ± 1.3) 
43.7 ±2.2(18.8 ±1.5' 
40.5 ±2.1(17.6 ±1.6 
9.7 ± 1.4(9.7 ± 1.4) 



9.9 ± 1.0(9.9 ± 1.0 
9.9 ±1.0(9.9 ±1.0 
10.0 ± 1.0(10.0 ± l.b) 
9.8 ±1.0(9.8 ±1.0) 

9.7 ± 1.0(9.7 ± I.O) 

9.8 ± 1.0(9.8 ±1.0) 

(TVq ~ lOkms^^ 
10.0'±Tl(10.0 ± 1.1) 

10.0 ±1.1(9.9 ±1.1) 

9.9 ±1.1(9.9 ±1.1) 

10.1 ± 1.1(10.1 ± 1.1) 

76.8 ±2.9(10.0 ±1.1) 

75.9 ±2.9(9.8 ±1.1) 
9.9 ±1.2(9.9 ±1.2) 



Y 


Y 


Y 


Y 


Y 


Y 


N 


Y 


N 


N 


N 


N 


Y 


Y 





= lOkms ^ 




10.1 ± 1.0 


10.1 ± 1.0 


Y 


Y 


10.2 ± 1.0 


10.2 ± 1.0 


Y 


Y 


10.2 ± 1.0 


10.1 ± 1.0 


Y 


Y 


10.1 ± 1.0 


10.1 ± 1.0 


Y 


Y 


11.0 ± 1.0 


11.0 ± 1.0 


Y 


Y 


10.5 ±1.0 


10.5 ±1.0 


Y 


Y 


10.0 ±1.0 


10.0 ± 1.0 


Y 


Y 




= lOkms ^ 




10.1 ± 1.1 


10.0 ±1.1) Y 


Y 


9.9 ± 1.0 


9.9 ± 1.0) 


Y 


Y 


9.9 ± 1.1 


9.9 ± l.l) 


Y 


Y 


10.2 ± 1.1 


10.2 ± 1.1 


Y 


Y 


37.6 ±2.0 


14.7 ± 1.3 


N 


N 


30.1 ± 1.8 


12.4 ±1.2 


N 


N 


9.9 ±1.2 


9.9 ±1.2) 


Y 


Y 





= lOkms ^ 






10.6 ± 1.4 


10.6 ± 1.4) 


Y 


Y 


11.3 ± 1.2 


10.3 ± l.l) 


N 


Y 


9.7 ± 1.2 


9.7 ± 1.2) 


Y 


Y 


32.8 ± 2.4 


11.5 ± 1.4) 


N 


N 


47.4 ± 2.3 


18.5 ± 1.6) 


N 


N 


46.3 ± 2.3 


16.4 ±1.7) 


N 


N 


10.3 ±1.5 


10.3 ± 1.5) 


Y 


Y 




= lOkms^^ 






10.0 ± 1.0 


10.0 ± 1.0) 


Y 


Y 



Y 


Y 


Y 


Y 


Y 


Y 


Y 


Y 


N 


Y 


N 


Y 


Y 


Y 



EM Algorithm 
TABLE 2 — Continued 



N 




^wrong 








Success? 








(km s-l) 




(km s ) 





Simulation N = 3000 


Nmem = 600 




{V)mem = 200kms-^ 


EM 


601(601) 


17(17) 


200.9 ±0.1(200.9 ±0.1) 


EMv 


597(596) 


76(76) 


200.6 ± 0.1(200.6 ±0.l) 


EMyfl 


597(597) 


59(59) 


200.8 ± 0.1(200.8 ±0.l) 


FMvw 


926(596) 


340(29) 


157.1 ±0.1(200.8 ±0.1) 


Za 


2982(722) 


2382(124) 


96.4 ±0.1(198.4 ±0.1) 


VT3,. 


2964(469) 


2364(167) 


96.6 ± 0.1(200.1 ±0.l) 




405(405) 


197(197) 


201.0 ±0.1(201.0 ±0.1) 





= lOkms"^ 




10.4 ±1.4 


10.4 ± 1.4 


Y 


Y 


10.2 ± 1.3 


10.2 ± 1.3 


Y 


Y 


10.2 ± 1.4 


10.2 ± 1.4 


Y 


Y 


70.8 ± 3.6 


10.4 ± 1.4 


N 


Y 


73.8 ±2.9 


13.7 ± 1.6 


N 


N 


73.9 ±2.9 


11.2 ± 1.6 


N 


Y 


10.3 ±1.6 


10.3 ± 1.6 


Y 


Y 



TABLE 3 

Results of Algorithms * with Simulated Data and o-Vq 



: 4 KM s 



N 




^ wrong 


{^) mem 








Success? 








(km S^^) 




(km s ) 







Simulation N = 30 

EM 
EMv 
EMvR 
EMvw 
3a 

VTsa 

VTgM 

Simulation N = 30 Nn 



24 



23 





0) 


25 


2 


2) 


24 


1(1) 


23 





0) 


25 


3 


3) 


23 


3 


3) 


23 


1 


1) 



15 



EM 


14(14) 


0(0) 


EMv 


20(14) 


12(2) 


EMvR 


15(14) 




EMvw 


14(14) 




3a 


30(21) 


\m 


VT3a 


26(18) 




VTbm 


13(13) 


2(2) 


Simulation N = 30 


Nmem = 6 




EM 


21(5) 


28(0) 


EMv 


6(3) 




EMvR 


108 


f, 


EMvw 


20(5) 


28(1) 


3a 


29(18) 


23(12) 




28(14) 


22(8) 




15(3) 


21(3) 



N = 30 Nn 



Simulation N = 30 N 

EM 
EMv 
EMvR 
EMvw 
3a 

VTsa 

VTbm 

Simulation 

EM 
EMv 
EMvH 
EMviv 
3a 

VTsa 

VTbm 

Simulation 

EM 
EMv 
EMvH 
EMviv 
3a 

VTsa 

VTbm 

Simulation 

EM 
EMv 
EMvH 
EMvw 
3a 

VTsa 

VTbm 

Simulation N = 30 N„ 



mem — 


24 




23 


23 








25 


25 


2 


2 


25 


25 


1 


1 


24 


24 


1 


1 


26 


26 


2 


2 


24 


24 


2 


2 


22 


22 


2 


2 



15 



N = 30 



15(14) 
14(14) 
13(13) 
15(14) 
29(20) 
28(15) 
13(12) 

Nmem — 6 

8(6) 
13(6) 
11(6) 
13(6) 
29(11) 
28(7) 
6(5) 



N = 30 N„ 



: 24 



1(0) 

1(1) 
1(1) 

1(0) 

14(5 

13(2 

2(3)- 



2(1) 


13 


2) 


13 


l' 


12 




23 


?! 


22 




0(1) 



23 








23 








23 








23 








23 


6 


1 


22 


3 


2 


22 


2 


2 



(V)„,er, 

50.2 ± 0.4i 
49.8±0.3i 

50.1 ±0.4i 
50.2±0.4i 
50.0±0.3i 
50.5±0.4i 
50.3±0.4( 

('V')mem - 

49.2 ± 0.4 

38.8 ± 0.5 

45.3 ±0.6 
49.1±0.4i 

48.5 ±0.5 

57.1 ±0.5 

49.3 ±0.5 

(V)mem - 

65.6 ±0.6 

36.9 ± 0.5 

50.7 ±0.6 

65.8 ±0.6 

60.2 ±0.5 

63.4 ±0.5 

70.3 ± 0.7i 

(V) mem — 
99.4±0.5i 
99.8 ± 0.4i 

99.5 ± 0.4i 
99.7±0.5i 
99.8 ± 0.4i 
99.3 ± 0.4i 
99.3±0.5i 



= 50kms~-^ 
50.2 ±0.4) 
49.8 ±0.3) 

50.1 ±0.4) 

50.2 ± 0.4) 
50.0 ± 0.3) 
50.5 ± 0.4) 

50.3 ± 0.4) 

= SOkms"^ 
'49.2 ±0.4) 
"48.3 ±0.4) 
'48.9 ± 0.4) 
'49.1 ±0.4) 
'50.6 ±0.5) 
'51.4 ±0.6) 
'49.3 ± 0.5) 

50kms~^ 
'48.9 ±0.8) 
'41.5 ±0.7) 
'48.4 ± 0.7) 
'49.0 ± 0.8) 
'50.0 ± 0.6) 
'52.3 ±0.7) 
'47.8 ± 0.7) 

100kms~^ 

99.4 ±0.5) 
'99.8 ± 0.4) 
'99.5 ± 0.4) 
'99.7 ±0.5) 
'99.8 ±0.4) 
'99.3 ± 0.4) 
'99.3 ±0.5) 



4kms 



(V)mem = lOOkms-1 
95.7 ±0.5(97.8 ±0.4) 

97.6 ± 0.4(97.6 ± 0.4) 
98.1 ±0.4(98.1 ±0.4) 

95.7 ±0.5(97.8 ±0.4) 
80.4 ±0.5(97.0 ±0.5) 
82.6 ±0.5(96.1 ±0.4) 
98.0 ± 0.4(97.8 ± 0.4) 



15 



(V) mem 
73.6 ± 1 

43.0 ± 0, 
36.6 ± 

43.4 ± 0, 

55.5 ±0 

53.1 ±0 
99.5 ±0, 

(V)mem — 

200.8 ± 0.4i 

200.8 ± 0.4i 
200.8 ± 0.4i 
200.8 ± 0.4i 

171.8 ± 0.5i 
186.0 ±0.5i 

200.9 ± 0.4i 

(V)mem — 



= lOOkms"^ 
4(99.3 ± 0.6) 
.9(98.4 ± 0.6) 
9(99.2 ±0.6) 
9(99.2 ± 0.6) 
7(107.2 ±0.9 
8(104.4 ± 1.1 
1.6(99.8 ± 0.6)' 

200kms~^ 

'200.8 ± 0.4' 
200.8 ± 0.4' 
200.8 ±0.4' 

200.8 ±0.4' 
'201.0 ±0.4' 

200.9 ± 0.4' 

;2oo.9 ± 0.4; 

200kms~^ 



3.2 ± 1.3 


3.2 ± 1.3 


Y 


Y) 


3.4 ± 1.2 


3.4 ± 1.2 


Y 


y) 


3.1 ± 1.2 


3.1 ± 1.2 


Y 


y) 


3.2 ± 1.3 


3.2 ±1.3 


Y 


Y) 


3.3 ± 1.1 


3.3 ± 1.1, 


Y 


y) 


2.9 ± 1.1 


2.9 ± 1.1 


Y 


y) 


3.3 ± 1.3 


3.3 ± 1.3* 


Y 


y) 



CTVq njem ~ 4kmS 

4.0 ± l"5(4.0 ± 1.5) 

40.5 ±4.7(4.1 ±1.5' 

13.8 ± 2.1(4.2 ± 1.6 
3.9 ± 1.5(3.9 ± 1.5) 

47.6 ±5.5(11.8 ±2.9 
39.0 ±5.8(11.2 ±3.0 

4.4 ± 1.7(4.4 ± 1.7) 

(TVo mem ~ 4kmS^"'^ 

53.0'±To(4.3 ± 2.5) 

33.7 ±5.1(11.0 ±3.5" 

17.0 ±4.4(13.1 ±3.6' 
53.3 ±8.1(4.3 ±2.5) 

48.1 ±7.3(20.4 ±4.8' 
45.7 ±7.3(18.5 ±5.2' 

36.9 ± 6.9(4.5 ± 2.6)' 



o-v, 
5.3±2.0i 
5.4± 1.9i 
5.4±2.0i 
5.5±2.0i 
5.5±2.0i 
5.2±2.0i 
5.5±2.1i 



= 4kms 
'5.3 ±2.0) 
'5.4 ± 1.9) 
'5.4 ±2.0) 
'5.5 ±2.0) 
'5.5 ±2.0) 
'5.2 ± 2.0) 
'5.5 ±2.1) 

= 4kms^-' 



^^0,mem 

8.5 ± 4.0(3.2 ± 1.3) 
3.4 ± 1.4(3.4 ± 1.4) 
2.9 ± 1.3(2.9 ± 1.3) 

8.6 ± 4.0(3.2 ± 1.3) 
39.4 ±6.8(10.5 ±3.5) 
38.3 ±6.9(6.7 ±2.1) 
3.2 ± 1.4(3.3 ± 1.4) 

'^Vo.mem = 4kmS-^ 

45.6 ±12.1(1.9 ±1.2) 

54.3 ±10.5(2.2 ±1.3) 

52.0 ±10.5(2.1 ±1.3) 

54.4 ±10.6(1.8 ±1.2) 

53.7 ±9.2(13.7 ±4.1 

53.1 ±9.2(10.5 ±3.5' 
1.9 ±1.3(2.0 ±1.4) 



<^Vo,„em = 4kms 
3.2 ± 1.5(3.2 ± 1.5) 

3.2 ± 1.5(3.2 ± 1.5) 
3.2 ± 1.5(3.2 ± 1.5) 
3.2 ± 1.5(3.2 ± 1.5) 
58.4 ±6.5(3.1 ±1.5 
42.0 ±6.4(3.1 ±1.5 
3.1 ±1.5(3.1 ±1.5) 
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TABLE 3 



Continued 





N 


^me77 


N wrong 


{V) 
(km 


mem 

s-i) 


(km s ) 


Success? 


EM 




14 


14 


0(0) 


201.1 ± 0.5 


201.1 ±0.5 


4.3 ± 1.8(4.3 ± 1.8) 


Y 


Y^ 




EMv 




14 


14 


0(0) 


201.1 ± 0.5 


201.1 ±0.5 


4.3 ± 1.8(4.3 ± 1.8) 


Y 


Y 




EMvR 




14 


14 


0(0) 


201.1 ± 0.5 


201.1 ±0.5 


4.3 ± 1.8(4.3 ± 1.8) 


Y 


Y 




EMvw 




13 


13^ 




200.6 ± 0.4 


200.6 ± 0.4 


4.0 ± 1.5(4.0 ± 1.5) 


Y 


Y^ 








30 


15; 




124.9 ± 0.5 


201.1 ±0.5 


81.4 ± 8.6(4.3 ± 1.8) 


N 


Y) 




VT3,. 




30 


12 




124.9 ± 0.5 


201.2 ±0.5 


81.4 ±8.6(4.7 ±2.0) 


N 


Y 




VTbm 




12 


12; 


3(3) 


201.2 ± 0.5 


201.2 ±0.5 


4.7 ± 2.0(4.7 ± 2.0) 


Y 


Y^ 





Simulation 

EM 

EMy 

BMvR 
EMyw 

3(T 

VT3,. 
VTem 



N = 30 Nn 



= 6 (V>mem = 200kmS-^ 

5(5) 0(0) 198.0 ±0.9(198.0 ±0.9 

5(5) 0(0) 198.0 ±0.9(198.0 ±0.9 

4(4) 1(1) 195.4 ±0.6(195.3 ±0.6 

5(5) 0(0) 198.2 ±0.9(198.2 ±0.9 

30(8) 24(2) 101.9 ±0.6(192.0 ±0.9 

30(5) 24(1) 101.9 ±0.6(195.2 ±0.6 

4(4) 2(2) 194.9 ±0.6(194.9 ±0.6' 



CTVo „eni ~ 4kmS 

3.6 ± 2"3(3.6 ± 2.3) Y(Y 

3.6 ±2.3(3.6 ±2.3) Y(Y 

0.4 ±0.6(0.2 ±0.4) N(N 

3.8 ±2.4(3.8 ±2.4) Y(Y 

65.4 ±9.1(11.8 ±4.8) N(N' 

65.4 ±9.1(0.1 ±0.2) N(N 

0.0 ±0.0(0.0 ±0.0) N(N 



APPENDIX 

A. ERRORS FOR THE ESTIMATES OF MEAN AND VARIANCE 

Suppose X is a random variable that follows a Gaussian distribution with mean {X) and variance cr^^. We sample X 
such that our data set of N observations is given by {{Xi, crjcJIi^i, where is the measurement error. We calculate 
estimates {X) and axo iteratively from Equations [TSlfTH To obtain the la errorbar associated with each estimate we 
propagate the measurement errors as well as the parameter errors from the previous iteration. The variance of any 
non-linear function of z variables, f{xi, X2, Xz), can be approximated by 



so long as xi...Xz are uncorrelated. From Equations fT3lfT4l the parameter estimates in iteration n - 



(Al) 



1 are functions of 



the data as well as the parameter estimates obtained the n*'* iteration. We must therefore propagate the variances cr^ 



n + 1 are 



and 



where 



and cr^{„}. From Equation lAU the variances associated with the parameter estimates obtained in iteration 
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^ 2 /^2{«}x3 ' 



f p{"} 
r - V 

,_i (tJ Ml + cr 



After 15-20 iterations, parameter estimates and the associated errors are insensitive to their (arbitrary) initial values 
(see Section S]). 

B. POOL-ADJACENT- VIOLATORS ALGORITHM FOR MONOTONIC REGRESSION 

Monotonic regression provides a nonparametric, least-squares fit to an ordered set of data points, subject to the 
constraint that the fit must be either non-increasing or non-decreasing. In our application, we have estimates 
of the probability that the i^^ star is a dSph member, and we wish to estimate a function p(a), the unconditional 
probability of membership as a function of elliptical radius a. We assume p(a) is a non-increasing function. If we sort 
data points by order of increasing a^, the expression for the monotonic regression estimate of p(a) is 

'I \ ■ ^j=uPM, 

piai) = mm max — , (ol) 

l<u<i i<v<N V — u + 1 

In the case that < Pau^i for all i, then we obtain the simple result p{ai) — Pa/. . However, when the data "violate" 

the shape restriction, i.e., when > Pmi^i for some i, monotonic regression determines the non-increasing function 
that provides a least-squares fit to the data. 

The "Pool- Adjacent- Violators" (PAV) algorithm calculates monotonic regression estimates simply and efficiently. For 
example, suppose we have a data set of = 6 stars and our ordered set of P^i is given by {1.0,0.9,0.8,0.5,0.6,0.2}. 
The adjacent probabilities {0.5, 0.6} "violate" our assumption that p{a) is non-increasing. The PAV algorithm identifies 
such "blocks" of violators and replaces them with the average value in the block. In this case, PAV would give 
p = {1.0, 0.9, 0.8, 0.55, 0.55, 0.2}. 

To implement monotonic regression we use the standard PAV algorithm of lGrotzinger fc Witzgalll (|1984f ). We search 
the ordered data points from left (smaller a.i) to right (larger a^) for violations between successive "blocks." When 
a violation is discovered, the offending points are "coalesced" into a single block having the average value of the 
contributing points. In the case that the newly formed block causes a violation with the block on its immediate 
left, these two blocks are coalesced into a still larger block, and so on until there is no violation with previously 
searched blocks. When no violation occurs, we advance the search again to the right. The resulting estimate, p(a), is 
a decreasing step function. Our PAV estimates for each galaxy are displayed in Figure [71 



C. RESULTS FROM EM ALGORITHM USING SIMULATED DATA 
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Fig. 9. — Performance of EM algorithm with simulated dSph data. Shown here are the artificial data and EM results for simulated samples 
of A'^ — 30 stars and assuming a dSph velocity distribution with mean 50 km s~^ and dispersion 10 km s~^. The three main panels depict 
simulations with different degrees of contamination: upper left, upper right, and lower left panels correspond to simulations with member fractions 
fmem — 0.2, 0.5, 0.8, respectively. Left sub-panels: angular distance (top) and magnesium index (bottom) versus velocity for the simulated data 
set. Text indicates the number of simulated data points, A'^, (including contamination), and the number, Nm, of those that are drawn from a 
dSph-like member population with mean velocity (V^)TTieTn and velocity dispersion (Tvq — 10 km s~^. Blue circles identify the simulated members; 
black squares identify the simulated contaminants. See Section[6]for further details of the dSph and contaminant distributions. Middle sub-panels: 
squares/circles represent the same simulated data, but color indicates the value of Pm resulting from the EM algorithm. Black (red; magenta; 
green; cyan; blue) markers signify Pm ^ 0.01 {Pm > 0.01; > 0.50; > 0.68; > 0.95; > 0.99). Dotted vertical lines enclose velocities that satisfy a 
conventional Sa clipping algorithm. In some of Figures [9]- 1201 one or both of these limits lies outside the plotting window. Text indicates estimates 
of the number of members, mean velocity and velocity dispersion returned by the EM algorithm. Right sub-panels: contaminants for which the 
EM algorithm gives Pm > 0.5 (black squares with X) or members for which it gives Pm < 0.5 (blue circles with X). 
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Fig. 10. — Same as Figure|9] but for simulated data sets with N — 30, {V)mem — 100 km s ^ and avQ,mem — 10 1cm s 
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Fig. 11. — Same as Figure|9] but for simulated data sets with N — 30, {V)mem — 200 km s ^ and avQ,mem — 10 1cm s 
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Fig. 12. — Same as Figure|9] but for simulated data sets with AT = 300, 



50 km s ^ and ctvq . 



10 1cm s~^ 
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Fig. 13. — Same as Figure|9] but for simulated data sets with N — 300, {V)mem — 100 l<m s ^ and avQ,mem — 10 km s 
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Fig. 14. — Same as Figure|9] but for simulated data sets with AT = 300, {V)r, 



200 l<m s ^ and ctvq, 



10 km s~^ 
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Fig. 15. — Same as Figure[9] but for simulated data sets with A'' — 3000, {V)mem — 50 km s ^ and (TVq, mem — 10 km s 
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Fig. 16. — Same as Figure|9] but for simulated data sets with N — 3000, {V)mem — 100 km s ^ and <7Vq, mem — 10 l^m s 
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Fig. 17. — Same as Figure|9] but for simulated data sets with N — 3000, {V)mem — 200 km s ^ and <7Vq, mem — 10 l^m s 
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Fig. 18. — Same as Figure|9] but for simulated data sets with N — 30, {V)mem — 50 1cm s ^ and avQ,mem — 4 km s 
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Fig. 19. — Same as Figure|9] but for simulated data sets with N — 30, {V)mem — 100 km s ^ and avQ,mem — 4 km s 
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Fig. 20. — Same as Figure|9] but for simulated data sets with N — 30, {V)mem — 200 km s ^ and avQ,mem — 4 km s 



